%# -*- coding: utf-8 -*-
%\documentclass[unicode,bachelor]{seuthesis} % 本科
\documentclass[master,oneside]{seuthesis} % 硕士
% \documentclass[doctor]{seuthesis} % 博士
% \documentclass[engineering]{seuthesis} % 工程硕士
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{amsthm}
 % 这里是导言区

\begin{document}
\categorynumber{000} % 分类采用《中国图书资料分类法》
\UDC{000}            %《国际十进分类法UDC》的类号
\secretlevel{公开}    %学位论文密级分为"公开"、"内部"、"秘密"和"机密"四种
\studentid{050962}   %学号要完整，前面的零不能省略。
\title{基于活动轮廓模型的超声图像分割}{}{Ultrasound Image Segmentation Using Active Contour Model}{}
\author{郭大威}{GUO Dawei}
\advisor{於文雪}{副教授}{YU Wenxue}{Associate Prof.}
%\coadvisor{副导师}{副教授}{Co-advisor's Name}{Associate Prof.} % 没有
                                % 可以不填
\degree{工学硕士} % 详细学位名称
\major[12em]{图像处理与科学可视化}
\defenddate{答辩日期}
\authorizedate{学位授予日期}
\department{计算机科学与工程学院}{School of Computer Science and Engineering}
\duration{2007.11—2008.6}
\address{群贤楼2楼}
\maketitle

\begin{abstract}{超声图像分割；活动轮廓模型；边界信息；区域信息；水平集方法}
\paragraph{论文题目：基于活动轮廓模型的超声图像分割方法研究}
\paragraph{学生姓名：郭大威}
\paragraph{导师姓名：於文雪}
\paragraph{学校名称：东南大学}

医学影像在医学辅助诊断和检测发挥着重要的作用，借助于为医生提供人眼所不能直接观察到的疾病信息，能使得医生对疾病做出准确的判断。超声成像以其无创、无痛和无放射等特点越来越成为临床广泛应用的首选方式。对超声图像的精确分割或定位是HIFU治疗中的非常重要的环节，是后期重建、配准以及治疗能达到良好效果的前提。本文针对这一主题进行了大量的研究工作，并将初步研究成果应用于实际的超声图像分割中，取得了一定的结果。

讨论了基于边界信息的水平集活动轮廓模型，包括几何活动轮廓模型和测地活动轮廓模型等。首先简要介绍了几何模型的优点和缺点，并通过分析Snake模型能量项推导出GAC模型，将能量函数与曲线演化相结合。引入了测地曲线计算的概念，解释了测地曲线长度和欧式曲线长度的联系和区别，跟据两个模型的最终演化方程，解释了GAC模型优于几何活动轮廓模型的原因，即边界吸引因子的存在。针对包括GAC模型在内的传统水平集方法需要进行重新初始化过程的问题，本文引入一种对GAC模型的改进模型，在模型中加入惩罚项和面积项（或加速项），从而完全消去了重新初始化的必要性，同时也加快了曲线的演化过程缩短了模型的演化时间。为了能保证结果的稳定性，讨论了该模型的数值实现以及时间步长的选择范围等重要问题。

讨论了基于区域信息的水平集活动轮廓模型，首先简要介绍了M-S分割模型理论，并通过该模型解释了M-S模型的一种简化模型，首个基于区域的活动轮廓模型C-V模型。该模型对初始化轮廓线不敏感，实际应用比较广泛。然而，该模型在不符合C-V模型的PC假设条件的情况下，很难得到良好的分割结果。针对C-V模型的PC条件假设前提的缺点，介绍了对该模型的一种改进模型RSF模型。RSF模型通过定义局部灰度拟合能量项，并且将局部的灰度变化作为分割算法的特征考虑，加入到总的能量方程中，所以该模型在不符合PC假设条件的图像分割中能得到较之C-V模型更好的分割结果。最后，本文提出了一种基于混合信息的活动轮廓模型图像分割方法。该方法在充分利用图像梯度信息的同时保留了区域信息来共同引导曲线的演化，因而相比单独依赖梯度或区域信息的方法在弱边界和层次丰富的超声图像能有较好的分割效果。
\end{abstract}

\begin{englishabstract}{Ultrasound Image Segmentation; Active Contour Model; Boundary Based; Region Based; Level Set Method}
    English abstract.
\end{englishabstract}

\begin{terminology}
  本论文专用术语的注释表
\end{terminology}

%----------------------------------------------开始正文-----------------------
\begin{Main} % 开始正文
%--------------------------------------------第一章-------------------------
\chapter{绪论}
医学影像向医生提供了人眼所不能直接观察到的人体内部结构和疾病信息，借助于这些信息医生可以更直观地检测和诊断病变的位置和形态等，进而对疾病做出准确的判断，因此，医学影像在临床上发挥着重要的作用。医学影像的成像技术（设备）利用物理学、电子学、计算机科学等基础科学的先进技术来诊断和治疗疾病，是一门多学科的交叉技术。自1895年伦琴发现X线不久，X线就被应用于人体疾病诊断检查，从此奠定了医学影像学的基础。从20世纪五六十年代开始又相继出现了超声成像、X线计算机体层成像（Computer Tomography，CT）、磁共振成像（Magnetic Resonance Imaging，MRI）和正电子发射体层成像（Position Emission Tomography，PET）等新的成像技术，并在医疗辅助诊断和疾病辅助治疗中充当重要的角色\cite{ZhuangTG1991}。医学影像学不仅扩大了人体的检查范围，提高了诊断水平，而且可以对某些疾病进行治疗\cite{WuEH2003}。

超声医学成像技术（设备）以超声波作为探测手段，通过接收人体脏器内部及表面的散射或反射波显示人体脏器的二维或三维图像。自20世纪五六十年代年商用超声产品问世以来，超声医学成像设备就以其无创、无痛、无放射线损伤、快速、准确、便捷、经济等特点成为临床诊断中重要的医学成像设备，是现代医学影像技术的主要组成部分同时也是临床最常用的检查方式之一。目前超声医学图像在腹部、心脏、皮肤浅表部位的疾病诊断以及胎儿的产前检查等领域得到了广泛的应用\cite{GaoSK2000}。

\section{超声图像成像原理}
超声是超过正常人耳能听到的声波，频率在20000赫兹（Hz）以上。超声检查是利用超声的物理特性和人体器官组织声学性质上的差异，以波形、曲线或图像的形式显示和记录，借以进行疾病诊断的检查方法。四十年代初就已探索利用超声检查人体，五十年代已研究、使用超声使器官构成超声层面图像，七十年代初又发展了实时超声技术，可观察心脏及胎儿活动。超声诊断由于设备不似CT或MRI设备那样昂贵，可获得器官的任意断面图像，还可观察运动器官的活动情况，成像快，诊断及时，无痛苦与危险，属于非损伤性检查，因此，在临床上应用已普及，是医学影像学中的重要组成部分\cite{WuEH2003}。不足之处在于图像的对比分辨力和空间分辨力不如CT和MRI高。

\subsection{超声的物理特性}
超声是机械波，由物体机械振动产生。具有波长、频率和传播速度等物理量。用于医学上的超声频率为2.5\~10MHz，常用的是2.5\~5MHz。超声需要在介质中传播，其速度因介质不同而异，在固体中最快，液体中次之，气体中最慢。在人体软组织中约为150m/s。介质有一定的声阻抗，声阻抗等于该介质密度与超声速度的乘积。

超声在介质中以直线传播,有良好的指向性.这是可以用超声对人体器官进行探测的基础。当超声传经两种声阻抗不同相邻介质的界面时其声阻抗差大于0.1\%,而界面又明显大于波长,即大界面时,则发生反射,一部分声能在界面后方的相邻介质中产生折射,超声继续传播,遇到另一个界面再产生反射,直至声能耗竭。反射回来的超声为回声。声阻抗差越大，则反射越强，如果界面比波长小，即小界面时，则发生散射。超声在介质中传播还发生衰减，即振幅与强度减小。衰减与介质的衰减系数成正比，与距离平方成反比，还与介质的吸收及散射有关。超声还有多普勒应（Doppler effect），活动的界面对声源作相对运动可改变反射回声的回率。这种效应使超声能探查心脏活动和胎儿活动以及血流状态。

\subsection{超声的成像基本原理}
人体结构对超声而言是一个复杂的介质，各种器官与组织，包括病理组织有它特定的声阻抗（\ref{tb:1-1}）和衰减特性。因而构成声阻抗上的差别和衰减上的差异。超声射入体内，由表面到深部，将经过不同声阻抗和不同衰减特性的器官与组织，从而产生不同的反射与衰减。这种不同的反射与衰减是构成超声图像的基础。将接收到的回声，根据回声强弱，用明暗不同的光点依次显示在影屏上，则可显出人体的断面超声图像，称这为声像图（Echogram）。

\begin{table}
    \caption{人体不同介质的声速与声阻抗}
    \label{tb:1-1}
    \begin{tabular}[t]{cll}
        \hline
        介质 &\vline 密度 &\vline 超声纵波速度（m/s）\\
        \hline
        空气 &\vline 0.001293 &\vline 332 \\
        \hline
        水   &\vline 0.9934   &\vline 1523 \\
        \hline
    \end{tabular}
\end{table}

人体器官表面有被膜包绕,被膜同其下方组织的声阻抗差大,形成良好界面反射, 声像图上出现完整而清晰的周边回声,从而显出器官的轮廓。根据周边回声能判断器官的形状与大小。超声经过不同正常器官或病变的内部，其内部回声可以是无回声、低回声或不同程度的强回声。

无回声：是超声经过的区域没有反射，成为无回声的暗区（黑影），可能由下述情况造成：1.)液性暗区：均质的液体，声阻抗无差别或差很小，不构成反射界面，形成液性暗区，如血液、胆汁、尿和羊水等。这样，血管、胆囊、膀胱和羊膜腔等即呈液性暗区。病理情况下，如胸腔积液、心包积液、腹水、脓液、肾盂积水以及含液体的囊性肿物及包虫囊肿等也呈液性暗区，成为良好透声区。在暗区下方常见回声增强，出现亮的光带（白影）；2.)衰减暗区：肿瘤，如巨块型癌，由于肿瘤对超声的吸收，造成明显衰减，而没有回声，出现衰减暗区；3.)实质暗区：均质的实质，声阻抗差别小，可出现无回声暗区。肾实质、脾等正常组织和肾癌及透明性变等病变组织可表现为实质暗区。

低回声：实质器官如肝，内部回声为分布均匀的点状回声，在发生急性炎症，出现渗出时，其声阻抗比正常组织小，透声性增高，而出现低回声区（灰影）。

强回声：可以是较强回声、强回声和极强回声。1.)较强回声：实质器官内组织致密或血管增多的肿瘤，声阻抗差别大，反射界面增多，使局部回声增强，呈密集的光点或光团（灰白影），如癌、肌瘤及血管瘤等；2.)强回声：介质内部结构致密，与邻近的软组织或液体有明显的声阻抗差，引起强反射。例如骨质、结石、钙化，可出现带状或块状强回声区（白影），由于透声性差，下方声能衰减，而出现无回声暗区，即声影（acoustic shadow）；3.)极强回声：含气器官如肺、充气的胃肠，因与邻近软组织的声阻抗差别极大，声能几乎全部被反射回来，不能透射，而出现极强的光带。

\section{图像分割的定义}
在医学图像诊断中，人们往往仅对图像中的某些部分感兴趣，比如对肿瘤的位置区域（位置信息）。这些部分常称为目标区域（其它部分称为背景），常对应于图像中特定的、具有独特性质的区域。这些特性可以是灰度、颜色、纹理等，目标可以对应单个区域，也可以对应多个区域。为了确定目标肿瘤的位置信息，需要将这些有关区域分离提取出来，在此基础上才有可能对目标进一步处理，如测量、分析和治疗。图像分割就是指把图像分成各具特性的区域并提取出感兴趣目标的技术。


借助于几何概念图像分割可以定义如下\cite{ZhangYJ2005}：
令集合R 代表整个图像区域，对R的分割可看作将R分成若干个满足以下条件的非空子集（子区域）$R_1, R_2, \ldots, R_n$，如下所示。

    a.)$\bigcup\limits^n_{i=1} R_i = R$；

    b.)对所有的i和j，$i\neq j$，有$R_i \cap R_j =\emptyset$；

    c.)对$i=1,2,\ldots, n$，有$P(R_i)=TRUE$；

    d.)对$i\neq j$，有$P(R_i\cup R_j)=FALSE$；

    e.)对$i=1,2,\ldots, n$，$R_i$是连通的区域.
其中，$(R_i)$是对所有在集合$R_i$中元素的逻辑谓词；$\emptyset$ 是空集。

上述条件a.)指出分割所得到的全部子区域的总和（并集）应能包括图像中所有的像素，或者说分割应将图像中的每个像素都分进某一个子区域中；条件b.)指出各个子区域是互不重叠的，或者说一个像素不能同时属于两个区域；条件c.)指出在分割后得到的属于同一个区域中的像素应该具有某些相同特性；条件d.)指出在分割后得到的属于不同区域中的像素应该具有一些不同的特性；条件e.)要求同一个子区域内的像素应当是连通的。

图像分割的文字定义可表述为：图像分割就是将给定的图像根据一定的规则划分成若干个互不相交的小区域的过程，小区域是某种意义下具有共同属性的像素的连通集合。这些规则可以依据图像的灰度、颜色、纹理等界定。

\section{常用的超声图像分割方法}

与CT、MR等其他医学图像相比，超声图像具有其自身的特点，如斑点噪声大、对比度低、存在与器官组织相关的纹理等。超声图像的这些特点由其成像原理所决定的：声像图以明暗黑白之间不相同的灰度来反映回声的有无和强弱，在无回声的组织上为暗区，强回声为亮区。超声图像所具有的特点也使得对超声图像的分割是一个有相当难度的问题。

目前，临床应用中的超声成像系统中使用比较多的的分割方式有基于专家手动分割的方法和基于阈值（Threshold）的方法。显然，专家手动分割方法的结果最令人满意，然而由于这样重复的工作以及时间上的耗费，使得该手工方法在实际应用中的工作效率特别低，受到很大的限制。简单的阈值分割法的算法实现很方便，并且快速、简单，然而由于超声图像中与器官相关的纹理信息以及图像存在的斑点噪声导致的图像灰度不均一使得该方法得不到良好的分割结果。

因此，在超声图像的自动或半自动的分割方法研究和临床应用就变得越来越重要，也正在成为超声图像分割的发展趋势。一直以来很多理论研究人员和临床工作者都在这方面做了大量的研究和应用工作，也取得了一定的研究成果。

常用的超声图像分割算法有：活动轮廓模型方法（AC），人工神经网络（ANN），马尔可夫随机场（MRF），离散动态轮廓模型（DDC）等。不同分割方法的所用的图像信息又包括：灰度分布、边界梯度、相位信息、纹理和形状信息等。接下来将根据不同的临床应用来分类介绍超声图像中所应依赖的分割方法。

\subsection{心脏病（Cardiology）}
超声心动图（Echocardiography）作为医学超声应用领域之一，主要用于诊断和评估缺血性心力衰竭等心脏疾病。与之相关的自动分割和跟踪左心室的文献也是很广泛的。常用的超声心动图分割方法有：活动轮廓模型\cite{mishra2003ga}，人工神经网络\cite{binder1999artificial}，基于灰度分布的统计模型\cite{mignotte2001endocardial}和Markov随机场\cite{friedland1989automatic}等。

Mishra等\cite{mishra2003ga}使用活动轮廓模型（AC）的方法解决心内膜边界探测（Endocardial border detection）问题，其优化过程使用遗传算法处理。该方法利用低通滤波和形态学方法用来初步估计初始化轮廓，使用非线性的灰度梯度映射建立能量方程，并通过不断的迭代得到最终的轮廓线。

ANN常被用于基于区域的分割中。例如，Binder等\cite{binder1999artificial}采用两层BackPropagation（BP）网络对样本集进行训练。首先，在单帧图像中进行径向搜索候选边界点，并使用时空（spatio-temporal）方法连接离散的轮廓点，进而得到心内膜边界。该方法用于舒张末期和收缩末期的SAX图像分割中，能够正确分割38副图像中的34副，并且分割结果与专家手工分割结果的相关系数为0.99。

\subsection{乳腺癌（Breast Cancer）}
乳房超声成像通常是作为诊断乳腺癌时乳房体检和乳房X光检查的辅助手段，与X线相比，超声检查乳腺的最大优点是能鉴别肿块是实性或囊性。超声不仅可以清晰地显示乳腺的解剖结构，还可以了解乳腺内有无肿块，为鉴别良性和恶性肿块提供了重要的信息\cite{ZhangXL2001}。常见的乳腺超声图像的分割方法有：神经网络\cite{chen2002diagnosis, huang2004watershed}，形变活动轮廓模型\cite{madabhushi2003combining, sahiner2004computerized}等。

以神经网络（NN）为基础的方法已被证明在乳腺癌超声图像分割中应用很广泛。其目标是根据一系列的输入特性，来做出分类的决策。Chen等\cite{chen2002diagnosis}提出了NN方法，其输入为：均方差比，自相关比，以及（Daubechies）小波系数的分布失真和多层神经网络感知器（其中一个隐层的网络通过反向传播误差训练）。通过242例数据（其中，161例良性肿瘤，81例癌症）检测，该方法达到了98.77\%的敏感度（sensitivity）和81.77\%的特异性（specificity）。作者强烈建议图像的纹理作为算法的重要组成部分，并使得该算法比较成功。

Huang和Chen\cite{huang2004watershed} 提出了一种集成神经网络分类的优点和一个分水岭分割特点的乳腺肿瘤超声图像轮廓检测方法。此方法的新颖之处在于其预处理过程，通过选择一个准确的标记该过程可以有效地帮助分水岭算法进行分割。为了能自动地从九个（组成）一组的预定义自适应滤波器，作者建议使用一个基于纹理的自组织映射（Self-Organizing Map，SOM）神经网络。

\subsection{前列腺（Prostate）}
经直肠超声成像（TRUS）一直是前列腺癌癌症诊断的一个必不可少的工具。前列腺的体积和边界信息在诊断、治疗和前列腺癌跟踪中发挥了重要作用。一般情况下，前列腺边界经常由横向平行的2-D切片沿前列腺的长度来描绘。这也导致前列腺边界检查的一些（半）自动方法的提出和发展。前列腺边界检测的方法有DDC模型\cite{ladak2000prostate}，活动轮廓模型\cite{yu2004segmentation}，统计形状模型\cite{shen2003segmentation}，形变模型\cite{singh1992image}等。

Ladak\cite{ladak2000prostate}提出了一种基于离散动态轮廓（DDC）的半自动2-D前列腺边界分割模型，该模型是在Lobregt和Viergever模型\cite{lobregt1995discrete}基础上发展而来，作者还将该模型扩展到了3-D模型\cite{wang2003semiautomatic}中。该方法需要手动选择四个（3-D为六个）前列腺边界上的特定位置的点，然后采用适当的插值算法完成初始模型（边界）。所选定的点在使用了方向信息的DDC形变过程中不变（基于前列腺组织比非前列腺组织超声图像颜色暗）。此方法也可以通过提供更多的预选边界点，然后再进行DDC形变算法。

Shen\cite{shen2003segmentation}提出了一种统计形状模型可以自动检测2-D前列腺边界。该模型通过在形状规则化步骤考虑探头位置和图像特征描述相对探头的不变性，将TRUS图像的自身性质考虑在内。在不同的方向和模量上Gabor滤波器被用于数据保真和仿射不变特性，用于在一个层次结构中捕捉前列腺的形状，被用在该模型中。该模型的主要贡献在于事先在TRUS图像中对前列腺的几何形状建模，通过考虑Gabor滤波器正交特性，会提高通过Gabor滤波器方式建模的前列腺边界。Shen\cite{zhan2006deformable}还提出了基于Zernike矩和SVM的边缘检测方法。

\subsection{妇产科（Obstetrics and Gynecology）}
超声在血管疾病诊断和妇产科检测中由于它的天然非介入特性使其优势明显且应用广泛。在妇产科超声图像分割提供了重要的测量数据，以评估胎儿的增长和胎儿畸形的诊断。由于在胎儿的脸部和周围的羊水之间的灰度对比明显，使得自动边界检测算法容易实现。多数自动分割算法都集中在卵泡测量和囊肿的检测和测量上。主要的分割方法有：分裂-合并法\cite{muzzolini1993multiresolution}，分水岭和阈值相结合法\cite{potocnik2002automated}和基于局部灰度熵\cite{zimmer1996distribution}的分割方法。

Zimmer\cite{zimmer1996distribution}等将局部灰度熵（local gray level entropy）作为常规变量，并在不同的假设局部灰度的PDF遵守Gaussian分布， Rayleigh分布， Weibull分布或Nakagami分布的实验。实验表明，局部灰度熵可应用于图像灰度直方图有多个不同的灰度峰值的区域分割中。提出了基于熵的阈值分割方法，即最小交叉熵阈值（Minimum Cross Entropy Thresholding）方法，其中，将分段变量由灰度和局部熵的线性组合取代。此二元分割算法应用于卵巢囊肿的超声影像的分割、定量分析和恶性卵巢群的检测中。

\section{本文的组织结构}
本文的研究重点是活动轮廓模型（Active Contour Model），通过对超声图像和图像分割进行比较全面的研究基础上，进行了以下几个方面的工作：
1.	 对传统的超声图像分割技术进行介绍，又总结了超声图像分割的新方法和最新的国内外研究现状，引出本文所用的基于活动轮廓模型的超声图像分割方法。分析了参数活动轮廓模型方法，总结了参数活动轮廓模型所具有的拓扑变化不独立性的弱点；并深入研究了几何活动轮廓模型方法。
2.	 针对传统的水平集活动轮廓模型方法在演化过程中需要进行重新初始化的问题，本文对测地活动轮廓模型进行推导，并通过增加一个新的约束项（惩罚项）解决该问题。并通过实验结果证明了超声图像分割质量的一致性。
3.	为了解决Chan-Vese模型方法对不符合该模型的分片同质（Piece Constant）假设条件的层次丰富图像的分割不理想的问题，本文引入一种使用轮廓线边界区域信息的改进模型。通过实验说明了该模型在Chan-Vese模型不能达到理想效果的图像中也能获得比较满意的分割结果。
4.	 最后，本文综合利用边界信息和区域信息提出了新的活动轮廓模型方法，即基于混合信息的活动轮廓模型方法。该方法不仅能保证良好的分割结果而且在不同的初始化条件下也能得到满意的分割结果，保证分割结果的鲁棒性。

基于上述的内容，本文章节结构安排如下：

第一章介绍了本文研究的背景与意义，概述了超声图像的成像原理和超声图像分割的研究现状及存在的问题。第二章介绍了活动轮廓模型和水平集方法的基础，包括参数的活动轮廓模型和水平集活动轮廓模型以及水平集的数值实现问题。

第三章研究了基于边界的活动轮廓（Active Contour）模型的分割方法。介绍了基于边界信息的测地活动轮廓模型，在此基础上引入对测地活动轮廓模型的改进方法，避免了该模型的重新初始化问题。

第四章研究了基于区域的活动模型的分割方法。介绍了基于区域灰度信息的Chan-Vese模型算法，并对此模型所存在的“同质区域”问题进行改进。综合三、四两章的研究基础，本文提出新的活动轮廓模型，该模型同时利用图像的边界和区域信息，在避免模型重新初始化的前提下，利用区域灰度信息，提高了超声图像分割的效果和鲁棒性。
%--------------------------------------------第二章-------------------------
\chapter{活动轮廓模型与水平集方法}
\section{活动轮廓模型概述}
自1987年Kass、Witkin和Terzopoulos首次提出活动轮廓模型，也称为Snake模型\cite{kass1988snakes}，并成功地应用于二维和三维图像分割问题，此后，活动轮廓模型方法逐渐成为图像分割领域最为活跃和成功的方法之一。Snake模型将一系列的图像处理（分割）问题转换为能量极小化问题。Snake模型利用图像的先验信息\cite{kass1988snakes}：在图像感兴趣的区域上给出一条封闭曲线作为初始轮廓曲线，该曲线通过能量函数控制，最小化能量函数使轮廓曲线在图像中发生形变，逐渐逼近待分割的目标边界。活动轮廓模型提出之前，MIT 人工智能实验室的Marr的图像处理分层理论认为，图像处理需要三个独立的层次来表达图像处理过程，人们只能依赖于从图像数据本身获得的底层信息，不可能使用各种先验知识等高层信息。这种严格的顺序研究方法将图像处理任务分成几个独立的阶段，但同时将底层的误差传播到了高层，没有修正的机会\cite{ZhouCX2005}。Snake模型能将待分割目标的知识和经验，如目标的形状，灰度、色彩等经验知识，和基于图像本身的低层视觉属性，如边缘、纹理、灰度等信息，有机地结合起来，得到待分割目标的完整表达。

自Snake模型提出以来，各种形式的活动轮廓模型如雨后春笋般地发展起来。根据轮廓线表示形式的不同，活动轮廓模型大致可以分为：参数活动轮廓模型和水平集活动轮廓模型\cite{li2005level}。参数活动轮廓模型直接以轮廓曲线的参数化形式表达轮廓曲线的形变过程，表达形式紧凑，有利于参数活动轮廓模型的快速实现，然而，参数活动轮廓模型的表达形式与轮廓曲线拓扑结构的变化不独立，比如曲线的合并或分裂等。基于曲线理论和水平集方法的水平集活动轮廓模型将轮廓曲线间接地表达为水平集函数的零水平集的形式，可以自然地处理轮廓曲线拓扑结构的变化。因而，水平集活动轮廓模型方法在实际的分割方法中得到广泛的应用。

与第一章所述的其它类分割算法相比较，基于活动轮廓模型（Active Contour-Based， AC Based）的图像分割算法具有以下几个特点\cite{ZhuGP2007}：

1.）活动轮廓模型的表达形式和实现形式多种多样。如上所述，活动轮廓模型按照轮廓线的表示方式可分为参数活动轮廓模型和水平集活动轮廓模型。参数活动轮廓模型中有经典活动轮廓（Snake）模型、气球（Balloon）模型和梯度矢量流（GVF）模型；水平集活动轮廓模型又包括几何活动轮廓（Geometric AC）模型、测地活动轮廓模型（Geodesic AC）以及 Chan-Vese（C-V）模型等。

2.）活动轮廓模型的能量方程的结构灵活多变。不同于Marr的图像处理分层理论，活动轮廓模型可以融合多种先验信息，如图像的边缘信息、图像的区域信息以及目标形状特征等。这种能融合多种先验信息的特征是其他几类图像分割算法所不具有的属性。而且，通常几种轮廓线表示形式相一致的活动轮廓模型通过简单的加性结合即可获得一种新的活动轮廓模型。

3.）适用范围相对其他分割方法较广泛。活动轮廓模型的能量方程的结构灵活多变的特性使得该类分割方法的适用性扩大。不仅在传统的医学图像分割，在分辨率低、斑点噪声大的超声图像上也能得到良好的分割效果。

基于活动轮廓模型的图像分割算法也有其自身的缺陷和不足，活动轮廓模型在图像分割前必须进行初始化。另外，参数活动轮廓模型不能自动处理轮廓线的拓扑变化，水平集活动轮廓模型在运算时间上耗费很多等。然而随着研究不断地深入，基于活动轮廓模型的图像分割算法也在不断地更新和完善。接下来初步介绍活动轮廓模型的几种形式。

\section{参数活动轮廓模型方法}
最早出现的一类活动轮廓模型就是参数活动轮廓模型（Parametric Active Contour），以经典的Snake模型为代表。Snake模型将一系列的图像处理问题统一转换为能量极小化问题，该模型是关于一条参数化闭合曲线的能量表达式，由来自演化曲线自身的内力和图像信息的外力构成，内力约束它的形状，外力引导它的行为，该参数化曲线通过极小化它所表示的能量函数来达到提取目标边界的目的。经典参数活动轮廓模型的局限性主要表现在模型外力的作用范围小、对轮廓初始位置敏感、不能够收敛到轮廓的凹陷区域、不能够处理拓扑变化、计算复杂度高等。针对Snake 模型的这些缺陷，研究人员对参数活动轮廓模型进行了外力场设计方面的研究并提出包括气球模型\cite{cohen1991active}和GVF模型\cite{xu1998snakes}在内的多种模型。此类方法已广泛应用于图像分割\cite{dubuisson1996vehicle}、运动跟踪\cite{LiXY2002}、三维重建等计算机视觉领域。

\subsection{经典参数活动轮廓模型（Snake模型）}
经典参数活动轮廓模型（Snake 模型）在图像中目标边界附近给出一条封闭的初始化轮廓曲线，$C(s)=(x(s),y(s)),s\in[0,1]$，轮廓曲线$C(s)$的能量由内部能量和外部能量两部分组成。极小化能量时导出内力场和外力场，轮廓曲线$C(s)$在自身的内力和由图像数据产生的外力共同作用下运动（发生形变），运动的最终结果是使内力和外力达到平衡状态，进而轮廓曲线与图像中目标边界相近或一致。Snake 模型定义参数活动轮廓模型的轮廓曲线$C(s)$的能量泛函为：
\begin{equation}
E_{Snake} = \alpha\int_0^1{{|C'(s)|}^2 ds} + \beta\int_0^1 {|C''(s)| ds} - \gamma\int_0^1 {g(|\nabla u(C(s))|)ds}
\end{equation}
式中的前两项分别用来表示曲线$C(s)$的弹性和刚性，其中，$\alpha$和$\beta$分别表示曲线$C(s)$的弹性系数和刚性系数，$\alpha$和$\beta$取与曲线参$s$无关的正常数。且常定义仅仅与曲线$C(s)$有关的能量为：
\begin{equation}
E_{Int} = \alpha\int_0^1{{|C'(s)|}^2 ds} + \beta\int_0^1 {|C''(s)| ds}
\end{equation}
由于$E_{Int}$不依赖其他外在的因素，常被称为内部能量函数，函数的内能保持了曲线的光滑性。相应地，第三项为外部能量函数，对于给定的图像，常定义为边界梯度算子，是一个依赖于图像特征的能量函数，因此，被称为曲线$C(s)$的外部能量函数：
\begin{equation}
E_{Ext} = \gamma\int_0^1{g(|\nabla u(C(s))| ds}
\end{equation}
函数$g$的定义为：
\begin{equation}
g(\nabla u(x,y)) = \frac{1}{1+{|\nabla G_\sigma(x,y)*u(x,y)|}^2}
\end{equation}
其中，$\nabla$表示空间梯度算子，$*$表示卷积算子，$G_\sigma$表示标准差为$\sigma$的二维高斯滤波器。从公式定义可以得到外部能量与图像灰度变化相关。在图像边界的灰度值变化较大，相应的图像梯度值也较大，外部能量则趋于变小；而在非边界区域灰度值的变化不大，所以图像的梯度的幅值很小，外能趋于变大使得曲线加快演化。

通过极小化Snake模型的能量函数$E_{Snake}$，得到对应曲线$C(s)$的偏微分方程为：
\begin{equation}
\frac{\partial X}{\partial t}=\frac{\partial}{\partial s}(\alpha\frac{\partial X}{\partial s})-\frac{\partial^2}{\partial s^2}(\beta\frac{\partial^2X}{\partial s^2})-\nabla E_{Ext}
\end{equation}
根据此演化方程，在目标边界的时候内力和外力达到平衡（方程右侧的数值为零）则停止演化。

Snake模型的缺点就在于高斯势能力只在边界周围很小的区域内有效，所以，只有在初始化轮廓线设置在图像边界临近区域内，才能达到理想的演化效果。

\subsection{气球模型（Balloon模型）}
Snake模型的高斯势能力的作用范围小且初始化轮廓线必须在边界附近的特点，使得Snake 模型的人工介入过多，因而制约了 Snake 模型在实际图像处理系统中的广泛应用。针对Snake模型的这一缺点Cohen 通过模拟气球的膨胀收缩原理于 1991 年提出了著名的气球（Balloon）模型\cite{cohen1991active}，该模型可以克服经典Snake模型的高斯势能力作用范围小的缺点。Cohen定义气球力为：
\begin{equation}
    F_B(C) = k_B \vec N(C)
\end{equation}
其中，$\vec N$为轮廓线$C$的内向单位法向量，$k_B$为气球力的常系数。当$k_B$为正数时，$F_B$驱动轮廓线向内收缩；而当$k_B$取负数时，气球力则驱动模型轮廓线向外膨胀。$F_B$是一个与图像特征无关的外力，轮廓线在图像域的任何区域都可以受到此力的作用，扩大了气球力的作用范围。对应的气球模型为：
\begin{equation}
\frac{\partial X}{\partial t}=\frac{\partial}{\partial s}(\alpha\frac{\partial X}{\partial s})-\frac{\partial^2}{\partial s^2}(\beta\frac{\partial^2X}{\partial s^2})-\nabla E_{Ext} + F_B
\end{equation}
其中，气球力的内向单位法向量$\vec N$通过计算切向量$\vec T$归一化后旋转90度得到；且轮廓线每一次变动气球力也都需要做相应的更新。切向量的计算公式为：
\begin{equation}
    \vec T(X_i^n)=\frac{X_{i+1}^n - X_i^n}{\parallel X_{i+1}^n - X_i^n\parallel}+\frac{X_i^n - X_{i-1}^n}{\parallel X_i^n - X_{i-1}^n\parallel}
\end{equation}
气球模型克服了Snake模型的缺点，气球力比高斯外力的作用范围大大提高了。但是气球力模型也有自身的缺点，首先气球力是一个单向的驱动力，即轮廓线只能膨胀或者收缩，所以轮廓线的初始化必须全部在带分割目标的内部或者外部；其次气球模型在初始化边界与目标边界交叉（重合）时会发生“泄漏（漏气）”现象。

\subsection{梯度矢量流模型（GVF模型）}
针对Balloon模型的单一方向和“漏气”现象的缺点，Xu和Prince等人于1998年提出了一种新的梯度矢量流（Gradient Vector Flow，GVF）模型\cite{xu1998snakes}。此模型基于一种新的静态外力，它不随时间变化，也不依赖于模型轮廓曲线的位置。GVF力是一种利用图像边缘导出的外力，但不同于Gaussian势能力，GVF将图像中目标边界的梯度映射到较远的范围，因而有更大的作用范围，不仅能作用在图像的边界还能作用于整个图像；相对于Balloon模型的优势在于，GVF是双向驱动的力。

GVF模型定义梯度矢量场$V(x,y)=[u(x,y),v(x,y)]$，定义能量函数并通过此函数求取$V$：
\begin{equation}
    E_{GVF}(V)=\int\int {\mu ({|\nabla u|}^2+{|\nabla v|}^2)+{|\nabla f|}^2({V-|\nabla f|}^2)} dxdy
\end{equation}
其中，第一项为矢量场$V$的平滑项，第二项为矢量场$V$的数据项，μ是一个正的权值，用来平衡上式积分号中第一项与第二项的作用。在图像的非边界区域，即$\nabla f$取值较小时，公式中的第一项决定积分值，将矢量场$V$在图像的边界处的分布向两边延拓产生一个缓慢变化的场；在图像的边界区域，即$\nabla f$取值较大时，公式中的第二项决定积分值，第二项起着主导作用，该项的作用是使$V$的值接近于$\nabla f$，使得积分项接近于边界梯度的效果。

根据该能量函数的Euler-Lagrange方程求解得到关于V的偏微分方程：
\begin{equation}
    \frac{\partial V}{\partial t}=\mu \nabla^2 V-(V-\nabla f){|\nabla f|}^2
\end{equation}
其中，$\nabla^2$表示Laplacian算子，$\nabla^2={\partial^2}/{\partial x^2}+{\partial^2}/{\partial y^2}$，根据此方程所得到的GVF场$V$，并用以代替Gaussian势能力。同Gaussian势能力类似，GVF场$V$也是静态力，在演化过程中不需要随着轮廓线的变化的更新。

\section{水平集方法}
水平集方法\cite{osher1988fronts}（Level Set）最初用于解决遵循热力学方程下的火苗的形状变化，是一种有效地实现曲线演化的方法。该方法将二维曲线嵌于三维空间曲面中，作为曲面的零水平集，从而曲线演化的过程变为对曲面的Hamilton-Jacobi方程的求解过程。水平集方法解决了曲线演化过程中的奇异性问题以及曲线的拓扑变化问题（即分裂或合并），近年来水平集方法逐渐成为研究的热点，并被广泛应用于图像分割中。

\subsection{曲线演化理论}
曲线演化理论是仅利用曲线的单位法向量和曲率等几何测度来研究曲线随时间形变的理论。在几何活动轮廓模型中，描述轮廓曲线几何特征的两个重要参数是单位法向矢量$\vec N$和曲率$\kappa$（而在参数活动轮廓模型中，轮廓线形变所依赖的曲线导数是一个与曲线参数有关的变量，不属于几何测度）。设一个闭合的动态曲线为$C(s, t), t>0$，该曲线的演化可以用下面的偏微分方程来描述：
\begin{equation}
    \frac{\partial C(p,t)}{\partial t}=\vec V = \alpha(p,t)\vec T + \beta(p,t)\vec N, \quad C(p,0)=C_0(p)
\end{equation}
其中，$\alpha$和$\beta$分别为切向和法向速率，即曲线C任意一点的速度$\vec V$分别沿曲线的切线方向和法线方向分解为两个速度。曲线的几何形状的变化只和法向速度分量$\beta$有关而与切向速度$\alpha$无关（$\alpha$仅影响曲线的参数化表示，不会改变其形状和几何属性），所以，曲线演化方程常表示为：
\begin{equation}
    \frac{\partial C}{\partial t}=\beta \vec N
\end{equation}

常见的曲线运动是曲率运动和常量运动，对应于两种最基本的曲线演化方式常量演化和曲率演化。
常量演化是指曲线运动的法向速度$\beta=c$，$c$为常系数，相应的演化方程为：
\begin{equation}
\frac{\partial C}{\partial t}=c \vec N
\end{equation}
常量演化可对应于形态学算子中的膨胀和腐蚀算子，$c>0$时曲线收缩，$c<0$时曲线扩张。且常量演化会导致曲线出现尖角和拓扑的变化。

曲率演化是指曲线运动的法向速度$\beta=a\kappa$，$a$为常系数，$\kappa$为曲率，相应的演化方程为：
\begin{equation}
\frac{\partial C}{\partial t}=a\kappa \vec N
\end{equation}
曲率演化使得任意简单闭合曲线在该偏微分方程作用下越来越简单，曲线的曲率的“振荡”幅度越来越小，且曲率的过零点越来越少，直到最后趋于圆形\cite{mokhtarian1992theory}。

\begin{figure}[htpb]
    \centering
    \includegraphics[scale=0.5]{figures/ch2/const.jpg}
    \caption{常量演化}
    \label{fig:2-1}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[scale=0.5]{figures/ch2/cur.jpg}
    \caption{曲率演化}
    \label{fig:2-2}
\end{figure}
从图（\ref{fig:2-1}）和图（\ref{fig:2-2}）两个图示可以看出，常量演化以均匀的速度进行膨胀和收缩运动，容易产生尖角；曲线进行曲率演化时，会变得越来越平滑并趋于圆形。因此，曲线的常值演化类似于仅在气球力作用下的参数活动轮廓模型轮廓线的运动过程；而曲线的曲率演化则类似于仅在内力作用下的参数活动轮廓模型轮廓线的运动过程\cite{ZhuGP2007}。

\subsection{水平集方法}
给定平面上的一条封闭轮廓曲线$C$，若$d(x,y)$是点$(x,y)$到轮廓曲线的垂直距离，称点$(x,y)$与距离$d(x,y)$之间的函数为符号距离函数（Signed Distance Function，SDF），记为$\phi(x,y)$。通常规定在轮廓曲线内部点$(x,y)$的符号距离函数$d(x,y)<0$为负值，在轮廓曲线外部点$(x,y)$的符号距离函数$d(x,y)>0$为正值，当$d(x,y)=0$的点的集合称为零水平集，它描述了平面上的一条封闭轮廓曲线。符号距离函数$\phi(x,y)$所取不同的数值的集合称为水平集（Level Set）。

水平集方法是由Osher和Sethian\cite{osher1988fronts}于1989年首次提出的一种非常有效的实现曲线演化的方法，将演化曲线$C(s ,t)$嵌入到一个比它高一维的函数$\phi(x ,y ,t)$中，其中函数$\phi$被称为水平集函数。通常水平集函数$\phi$用其零水平集来表示曲线$C$，即$C(s ,t)=\{(x ,y)|\phi(x ,y ,t)=0\}$。传统水平集函数一般被定义为如下符号距离函数（SDF）：
\begin{equation}\label{xx}
    \phi(x,y,t)=\left\{
        \begin{aligned}
            -d((x,y),C(s,t)), & \quad \text{如果点(x,y)在曲线C(s,t)内部;}\\
            0, & \quad \text{如果点(x,y)在曲线C(s,t)上;}\\
            d((x,y),C(s,t)), & \quad \text{如果点(x,y)在曲线C(s,t)外部.}
        \end{aligned}
    \right.
\end{equation}
其中$d((x, y), C(s, t))$表示坐标$(x, y)$到曲线$C(s, t)$的距离。

水平集方法将曲线的演化问题转化成了高维函数的更新问题。当高维函数$\phi$发生变化时，嵌入其中的曲线$C$也随之变化；一旦高维函数$\phi$被确定，则曲线$C$也最终被确定。假设曲线$C(s ,t)$的演化方程被定义为：
\begin{equation}
    \frac{\partial C}{\partial t}=V \vec N
\end{equation}
其中$V$表示曲线的演化速度，$\vec N$表示向内单位法向量。正如上一节所述，沿曲线切线方向的速度并不影响曲线在变形中的几何特征，因此，由公式所定义的演化方程的形式具有一般性。给定一个水平集函数$\phi(x, y, t)$，用其零水平集来表示演化曲线$C(s ,t)$，则有：
\begin{equation}
    \phi(C(s,t),t) = 0
\end{equation}

此轮廓线需满足曲线运动的偏微分方程，即${\partial C}/{\partial t}=V(C)\vec N$，对公式求全微分得到方程：

\begin{equation}
\frac{\partial \phi}{\partial t} + \nabla \frac{\partial C}{\partial t}=0
\end{equation}
其中，$\nabla \phi$是$\phi$的梯度，该公式为轮廓线对应的水平集的演化方程。

将当前的轮廓曲线$C$嵌入到一个更高维的水平集函数$\phi$中，并利用轮廓曲线运动方程与 Hamilton-Jacobi 方程的相似性，水平集方法给出了一种轮廓曲线运动的强鲁棒性的计算方法。水平集方法最大的优势在于它的稳定性以及拓扑无关性。

\section{水平集活动轮廓模型方法}
水平集方法和活动轮廓模型相结合的曲线演化方法是目前应用广泛的一种图像分割方法。水平集活动轮廓模型方法利用曲线的几何特征，通过建立与轮廓线相关的能量函数方程，最小化该能量函数方程，并使得曲线进行演化（发生形变）逼近图像中的目标边界。曲线的演化方程最终转化为水平集函数的偏微分方程问题。水平集活动轮廓模型又可以分为基于边界\cite{caselles1995geometric,li2005level,caselles1997geodesic}和基于区域\cite{ chan2001active, wang2010efficient}等几类活动轮廓模型。

\subsection{几何活动轮廓模型}
以 Snake 模型为代表的参数活动轮廓模型不能自动地处理轮廓线的拓扑变化。水平集活动轮廓模型利用水平集函数来隐含地表示模型轮廓线，从而成功地解决了轮廓线拓扑变化的问题。几何活动轮廓模型（Geometric Active Contour Model）是由 Caselles 等人与 Malladi\cite{caselles1995geometric, malladi1995shape}等人分别独立地提出。几何活动轮廓模型的提出开创了一种新形式的活动轮廓模型——水平集活动轮廓模型。

对于一条运动轮廓曲线$C(s,t)=(x(s,t),y(s,t))$，其中$s$表示曲线的参变量，$t$表示时间变量，$C$表示几何活动轮廓模型的轮廓线，则几何活动轮廓模型的轮廓线$C$的演化方程可表示为：
\begin{equation}
    \frac{\partial C}{\partial t} = g(V+\kappa)\vec N
\end{equation}
其中$\vec N$为曲线$C$的内向单位法向量，$\kappa$为曲线$C$的曲率，$V$为常系数（常量演化速度），$g$为依赖于图像特征的函数（即边界指示因子）被定义为：
\begin{equation}
    g=\frac{1}{1+{(G_{\sigma}*u)}^2}
\end{equation}
其中，$\nabla$表示空间梯度算子，$*$表示卷积算子，$G_{\sigma}$表示标准差为$\sigma$的二维高斯滤波器。边界指示因子的作用可解释为：当坐标点$(x,y)$在图像边界周围时，图像梯度的值就相对很大，所以$g(x,y)$的值会约等于$0$；相反，当坐标点$(x,y)$远离图像边界时，图像梯度的值就相对很小，所以$g(x,y)$的值会约等于$1$。由此，可以通过$g(x,y)$的值来确定对应像素点是否位于图像的边界区域，当轮廓线演化到图像边界时则会停止演化。通过曲线演化方程和水平集方法，几何活动轮廓模型的曲线演化方程可转为如下形式：
\begin{equation}
    \frac{\partial \phi}{\partial t}=g|\nabla \phi|(V_0+div(\frac{\nabla \phi}{|\nabla \phi|}))
\end{equation}
其中，$\phi$为水平集函数，对应的零水平集为轮廓曲线$C$，通过此水平集的演化最终可以达到分割的效果。

几何活动轮廓模型的设计直接基于曲线演化理论和水平集函数方法，没有将类似传统的参数活动轮廓模型的能量综合，因此很难成为一种有效的水平集模型方法。
\subsection{测地活动轮廓模型}
本节将介绍另一种经典的水平集活动轮廓模型——测地活动轮廓模型（Geodesic Active Contour Model）\cite{caselles1997geodesic}。测地活动轮廓模型由原始的Snake模型推导而来，即通过设计轮廓线的能量泛函，极小化能量泛函得到轮廓线的运动方程。与几何活动轮廓模型相比较，设计过程比较复杂，但为研究者们提供了一种更开阔的设计思路\cite{ZhuGP2007}。下一章将详细推导和解释测地活动轮廓模型。

令演化曲线$C(p, t)$表示测地活动轮廓模型的轮廓线，其中$p\in [0,1]$表示曲线$C$的参变量，$t$表示时间变量。在Riemannian空间新的距离定义为：
\begin{equation}
    L_R=\int_0^L{g(|\nabla I(C(q))|)|C'(q)|}dq=\int_0^L{g(|\nabla I(C(q))|)}ds
\end{equation}
该公式定义了Riemannian空间的曲线长度。对比欧式空间的曲线长度，该定义增加了一个指示图像边界的权重系数$g(|\nabla I(C(q))|)$。根据最陡（梯度）下降法和$L_R$的Euler-Lagrange方程可得到该能量函数的演化方程为：
\begin{equation}
    \frac{\partial C(t)}{\partial t}=g(I)k\vec N-(\nabla g\cdot \vec N)\vec N
\end{equation}
其中，$k$为曲线的欧式曲率，$\vec N$为单位内向法向量公式提供了一个最小化$L_R$曲线的曲线演化流。根据曲线演化方程和水平集方法，得出以下水平集方程：
\begin{equation}
    \frac{\partial u}{\partial t}=|\nabla u|div(g(I)\frac{\nabla u}{|\nabla u|})=g|\nabla u|k+\nabla g \cdot \nabla u
\end{equation}

相比于几何活动轮廓模型的水平集演化方程，测地活动轮廓模型的水平集演化方程多了最后一项$\nabla g \cdot \nabla u$，通常称之为吸引因子。当嵌入在$u$中的轮廓线$C$远离图像边界时或者与边界发生交叉时，该项产生的力把轮廓线$C$拉回边界的中间。因此，与几何活动轮廓相比而言，测地轮廓模型能防止轮廓线出现边界泄漏现象。

\section{本章小结}
本章介绍了参数活动轮廓模型的基本原理和三种经典的模型，包括Snake模型，气球模型和GVF模型，每一种模型都是对前一种模型的改进和完善。同时介绍了曲线演化理论和水平集方法，并且基于此介绍了水平集活动轮廓模型，包括几何活动轮廓模型和测地轮廓模型。水平集活动轮廓模型的优点在于其克服了参数活动轮廓模型的拓扑变换不独立的缺点。

%--------------------------------------------第三章-------------------------
\chapter{基于边界信息活动轮廓模型的超声图像分割}
根据活动轮廓模型的表示方式和实现方式，可以将其广义地分为两类：参数活动轮廓模型和水平集活动轮廓模型。参数活动轮廓模型显式地表示为在Lagrangian框架中的参数化曲线，而水平集活动轮廓模型被隐式地表示为在Eulerian框架中的演化的二维函数的水平集\cite{li2005level}。水平集活动轮廓模型又可以分为基于边界和基于区域等几类活动轮廓模型。本章将详细介绍基于边界信息的水平集活动轮廓模型，及该类方法对超声图像分割中的应用，下一章将详细讨论基于区域信息的水平集活动轮廓模型及其在超声分割中的应用。

本章结构安排如下：第一节简要介绍第一个水平集活动轮廓模型，几何活动轮廓模型，以及该模型的设计模式；第二节由Snake模型为基础推导出测地活动轮廓模型。第三节针对测地活动轮廓模型需要重新初始化过程的问题，给出对测地活动轮廓模型的一种改进方法，即无需重新初始化的活动轮廓模型。

\section{几何活动轮廓模型的提出}
几何活动轮廓模型分别由Caselles\cite{caselles1995geometric}和Malladi\cite{malladi1995shape}等独立引入。该模型根据平面曲线演化理论和水平集方法，将变化的轮廓线嵌入到更高维的隐函数中，通过对高维隐函数的演化相应的零水平集曲线也得到演化。

曲线演化可分为基于曲率的和基于常量的演化方法。基于曲率的演化方向和速度由曲线的法向和平均曲率决定。几何活动轮廓模型是基于以下的曲率演化方程：
\begin{equation}
    \frac{\partial \phi}{\partial t}=|\nabla \phi|(div(\frac{\nabla \phi}{|\nabla \phi|}))
\end{equation}

向该演化模型中添加常量演化项，此项表示一个指向水平集法向方向的恒定力；另外，为了能产生在目标边界停止的效果，将其与函数g(x)相乘得到：
\begin{equation}
\left\{ \begin{aligned}
    &\frac{\partial \phi}{\partial t}=g|\nabla \phi|(v+div(\frac{\nabla \phi}{|\nabla \phi|}))\\
    &\phi(0,x)=\phi_0 (x)
\end{aligned} \right.
\end{equation}
定义$g$为：
\begin{equation}
    g=\frac{1}{1+{(\nabla G_{\sigma}*u)}^2}
\end{equation}
其中，$v$是一个正常量，$\kappa=div(\nabla \phi/|\nabla \phi|)$是曲面的曲率，$\nabla G*u$是图像$u$与标准方差为$\sigma$的高斯函数$G_\sigma$的卷积，$u_0$为原始数据。当轮廓线在真实边界线附近即图像边缘时，$g(x)$会减小而趋于$0$进而使得曲线停止演化；相反，当轮廓线远离图像边缘时，$g(x)$会增大趋于$1$进而根据法向方向进行演化，因此，能够达到在目标边界停止的效果。

对应于该模型的轮廓线$C$的演化方程为：
\begin{equation}
    \frac{\partial C}{\partial t}=g(v+\kappa)\vec N
\end{equation}
其中，$\vec N$表示曲线$C$的内向法向向量，$\kappa$为曲线$C$的曲率。从该方程可以看出，曲线以$g(v+\kappa)$的速度进行演化：当$v>0$时轮廓线$C$向外进行膨胀演化；反之，当$v<0$时轮廓线向内进行收缩演化。

该模型的意义在于其开拓了一种新的活动轮廓模型，即水平集活动轮廓模型。该模型直接基于曲线演化方程设计（也成为轮廓演化模型），与参数活动轮廓模型相比，模型对拓扑结构的变化是独立的。

\section{测地活动轮廓模型的原理与实现}
测地活动轮廓模型不同于几何活动轮廓模型，该模型是将基于能量的活动轮廓（Snake）与测地计算（Riemmanian空间的最小化距离曲线）相关系的活动轮廓模型。从这个测地模型中推导出来活动轮廓模型的测地偏微分方程（Geodesic PDE）与之前的几何活动轮廓模型（轮廓演化模型）有很大的提高。

\subsection{基于能量的活动轮廓模型}
令$C(q):[0,1]$表示一个参数化的极坐标曲线，$I:[0,a]*[0,b]$为待检测边界的给定图像。与曲线$C$相关的传统的Snake模型方法如下：
\begin{equation}\label{eq:3-1}
    E(C)=\int_0^1{(\alpha|C'|^2+\beta|C''|^2 )}dq-\lambda\int_0^1{|\nabla I(C)|}dq
\end{equation}
其中，前两项控制待检测目标边界的光滑性（内能），第三项用于将轮廓线吸引到目标边界（外能）。解决Snake模型问题相当于，在给定$\alpha, \beta \text{和}\lambda$条件下，寻找一个能最小化能量$E$的曲线$C$。当初始化曲线$C$包围多个待分割物体，而检测到（分割）这些物体的可能性将不大，在此情况下就必须加入特殊的拓扑处理过程（Topology-Handling Procedures）。注意到这样的方式不是拓扑独立的，当该方法处理包含多个待分割物体的图像时将会比较麻烦。

Snake模型还可能存在的一个问题是需要选择参数来平衡边界光滑性和边界的准确性。考虑当刚性系数为零的时候，即当$\beta=0$时的一种特种情况下的Snake模型。这样假设的原因是：首先，不存在高阶导数的情况下，可以很自然地保持轮廓曲线的光滑性；其次，这样的参数设置可以推导出基于能量的活动轮廓模型和几何曲线演化直接的关系。当$\beta=0$时上述能量方程的第一项趋向于消除了刚性的存在。另外，规则化的GAC模型是从仅有第一项的模型中推导出来的。

设置参数$\beta=0$，同时替换边界检测算子$|\nabla I|$为更广义的梯度函数$g(|\nabla I|^2)$，即当$r\rightarrow\infty$时，$g(r)\rightarrow 0$的单调递减函数，此时得到如下模型：
\begin{equation}\label{eq:3-2}
    E(C)=\alpha\int_0^1{|C'|^2}dq+\lambda\int_0^1{g(|\nabla I(C)|^2)}dq
\end{equation}

注意到，最小化上述能量，相当于将曲线演化到$|\nabla I(C)|$最大的轮廓点上（边界检测算子），同时也保持曲线的光滑性。这样的效果和Snake模型中的效果是相同的。边界光滑性和边界的准确性之间的平衡由上述参数中的自由参数控制。文献\cite{caselles1997geodesic}推导并证明了，极小化上述能量模型可转化为如下形式：
\begin{equation}\label{eq:3-3}
    Min\int_0^1{g(|\nabla I(C(q))|)|C'(q)|}dq
\end{equation}
此时，已经将最小化公式（\ref{eq:3-2}）问题转化为新测量标准的Riemannian空间的测地曲线计算。

\subsection{测地曲线计算}
Riemannian空间和欧式空间的曲线长度对比可通过图示说明如下：
\begin{figure}[htpb]
    \centering
    \includegraphics[scale=0.5]{figures/ch3/geo.jpg}
    \caption{测地曲线与欧式曲线}
\end{figure}

不同的距离定义空间的最短距离是不同的，$E_{AB}$为欧式空间的最短距离，$G_{AB}$为Riemannian空间的最短距离。Riemannian空间的最短距离即定义为测地线，最小化Riemannian距离的过程即为测地线计算。

测地曲线的原理可做以下推导和解释。曲线$C$的欧式长度可以通过以下公式计算得到：
\begin{equation}\label{eq:3-4}
    L=\oint{|C'(q)|}dq=\oint ds
\end{equation}
其中，$ds$是欧式空间的弧长。很明显沿着矢量流$C=\kappa \vec N$（$\kappa$为欧式曲率）演化，是缩短$L$长度的最快方式，即沿着方程$L$的梯度方向演化。该矢量流常被称为欧拉曲线最短下降流（Euclidean Curve Shortening Flow）。而公式（\ref{eq:3-3}）是一种在Riemannian空间的新的曲线长度定义为：
\begin{equation}\label{eq:3-5}
    L_R=\int_0^L{g(|\nabla I(C(q))|)|C'(q)|}dq=\int_0^L{g(|\nabla I(C(q))|)}ds
\end{equation}
对比欧式空间的曲线长度公式，该定义在欧式长度微分$ds$上增加了一个指示图像边界的权重系数$g(|\nabla I(C(q))|)$。因此，此方法在探测物体时不仅仅能找到传统的最短路径，而且能够将图像的特征信息（比如边界）作为评判的标准。另外，注意到对函数$g$的假设仅限于递减的正函数，因此，基于测地计算的边界检测方法可以使用任意边界检测算子函数$g$。
\subsection{GAC模型的水平集演化}
为了最小化$L_R$，通过最陡（梯度）下降法找到$L_R$的梯度下降方向，因此，需要计算Euler-Lagrange方程。最陡（梯度）下降法，将初始化曲线$C(0)=C_0$演化到$L_R$的局部最小值，可以通过以下曲线演化方程：
\begin{equation}\label{eq:3-6}
    \frac{\partial C(t)}{\partial t}=g(I) \kappa \vec N - (\nabla g \cdot \vec N) \vec N
\end{equation}
其中，$\kappa$是欧式曲率，$\vec N$是内向单位向量，等号右边部分为公式（\ref{eq:3-3}）的Euler-Lagrange方程。该方程表明为了缩短$L_R$的长度，曲线$C$上的每一个点的运动方向和速度等。待分割物体就是在公式（\ref{eq:3-6}）稳定状态时候的结果。



根据曲线演化理论和水平集方法，测地线计算问题可以等价于$u$的（最陡）梯度下降演化，关于$u$的水平集方程为：
\begin{equation}\label{eq:3-7}
    \frac{\partial u}{\partial t}=|\nabla u|div(g(I)\cdot{\frac{\nabla u}{|\nabla u|}})=g|\nabla u|\kappa+\nabla g\cdot\nabla u
\end{equation}

其中，$u$是将轮廓线$C$嵌入的水平集，等式右边是关于$u$的Euler-Lagrange方程，$\kappa$通过水平集$u$计算得到，$\kappa=div(\nabla u /|\nabla u|)$。

在GAC模型中，曲线通过新的梯度项$\nabla g$产生的力被吸引到边界的中间（Middle of the Boundaries），$\nabla g \cdot \nabla u$称为吸引因子。当嵌入在$u$中的轮廓线$C$远离图像边界时或者与边界发生交叉时，该项产生的力把轮廓线$C$拉回边界的中间。因此，与几何活动轮廓相比而言，测地轮廓模型能防止轮廓线出现边界“泄漏”现象。

\section{无需初始化的活动轮廓模型}
对于给定的水平集函数，其演化方程可以表示为：
\begin{equation}\label{eq:3-8}
    \frac{\partial \phi}{\partial t} + F|\nabla \phi|=0
\end{equation}
其中，$F$为速度函数，此方程通常被称为水平集方程（Level Set Equation）。对于图像分割来说，速度函数$F$依赖于原始图像和水平集函数$\phi$。在传统的水平集方法中，水平集函数在演化的过程中可能产生尖点和（或）角点，这样会导致接下来的计算有很大的误差。

任何一个符号距离函数（Signed Distance Function，SDF）都必须满足一个特性，即$|\nabla \phi|=1$；作为SDF的水平集函数来讲也必须满足此条件。为了防止这些问题的出现，通常会在计算过程中加入一种数值计算策略（Scheme）：在函数演化之前，将其初始化为一个符号距离函数（SDF），然后，在演化过程中，将函数$\phi$再次“重新初始化（Re-initialize）”为一个符号距离函数（SDF）。重新初始化过程在传统的水平集方法中是很重要，而且是必须的。

\subsection{重新初始化的缺点}
在传统的水平集方法中，重新初始化作为一种数值弥补方法被广泛应用。标准的重新初始化方式通过解答以下的重新初始化方程：
\begin{equation}\label{eq:3-9}
    \frac{\partial \phi}{\partial t}=sign(\phi_0)(1-|\nabla \phi|)
\end{equation}
其中，$\phi_0$是需要重新初始化的函数，$sign(\phi)$是符号函数。已有很多关于重新初始化研究的文献资料\cite{peng1999pde}，而且其中的大部分都是基于PDE方法的变型。但是，当$\phi_0$不光滑或者$\phi_0$在曲面的一侧很陡峭（Steeper）在曲面的另一面却相对光滑，水平集函数$\phi$的零水平集的变换将会有偏差。而且，当水平集函数与SDF（特性）相差很远时，那些重新初始化方法也不足以将水平集函数重新初始化为SDF。在实际应用中，经过几次迭代后，演化的水平集函数将会偏离SDF（特性），尤其在时间间隔设置为很大的时候。

目前，重新初始化已被广泛用来作为维持曲线演化的稳定性和确保理想的演化结果数值补救策略（Numerical Scheme）。从实际应用的观点而言，重新初始化过程可能相当复杂，耗时，而且有些副作用。此外，大部分水平集方法也有自身的缺点，例如，何时以及如何重新初始化水平集函数为一个SDF\cite{gomes2002reconciling}。重新初始化过程很重要而且在基于PDE的方法中也是必须的，但是可以通过有限差分法消除。

\subsection{包含惩罚项的变分水平集模型}
任何一个符号距离函数（SDF）都必须满足一个特性，即$|\nabla \phi|=1$；相反，满足此特性的函数$\phi$都是一个SDF（加上一个常数）\cite{li2005level}。显然，在二维空间$\Omega \in \Re ^2$中，可以通过以下积分来测量一个函数$\phi$符合SDF特性的程度：
\begin{equation}\label{eq:3-10}
    P(\phi)=\int {\frac{1}{2} {(|\nabla \phi|-1}^2)} dxdy
\end{equation}
定义上述的能量项为惩罚项，综合此惩罚项建立变分方程如下：
\begin{equation}\label{eq:3-11}
    \xi(\phi)=\mu P(\phi)+\xi_m(\phi)
\end{equation}
其中，$\mu>0$为控制惩罚项（使得$\phi$趋近于SDF）系数，$\xi_m$为驱动函数$\phi$的零水平集的能量项。为了将此变分方程应用与图像分割中：能量项$\xi_m$被定义为依赖于图像数据的能量项（作为外能，external energy）；相应地，能量项$P$为$\phi$的函数且仅仅依赖$\phi$函数本身（作为内能，internal energy）。

在图像分割过程中，主动轮廓模型的轮廓曲线向目标边界移动。为实现这种演化，明确定义外部能量，使得零水平集曲线可以向目标边界移动。令$I$为原始图像，定义$g$为边界指示因子：
\begin{equation}\label{eq:3-12}
    g=\frac{1}{1+{(\nabla G_{\sigma}*I)}^2}
\end{equation}
其中，$G$是标准差为$\sigma$的高斯核函数。定义有关函数$\phi$的外部能量为：
\begin{equation}\label{eq:3-13}
    \begin{aligned}
        \xi_m(\phi)&=\lambda L_g(\phi)+\nu A_g(\phi)\\
        L_g(\phi)&=\int {g\delta(\phi)|\nabla \phi|}dxdy\\
        A_g(\phi)&=\int {gH(-\phi)}dxdy
    \end{aligned}
\end{equation}
其中，$\lambda>0$，v为常数，$\delta$是单变量的狄拉克函数（Dirac Function），$H$是传递函数（Heaviside Function）。

根据上节的GAC模型的解释，不难发现$L_g(\phi)$为计算函数$\phi$的零水平集曲线的长度。令$\phi$的零水平集可通过参数化表示为$C(p),p\in[0,1]$，则有$ds=g(C(q))|C'(q)|dq$；而$A_g(\phi)$的引入则是为了加快曲线的演化。注意到，当函数$g$是常数$1$时，$A_g(\phi)$为零水平集内部区域的面积；当函数$g$不为常数时，A$_g(\phi)$可以被作为是加权的零水平集内部区域面积。

系数$\nu$可以是正数或者负数，取决于初始轮廓线相对于目标对象的位置。例如，如果最初轮廓在目标对象的外围，加权区域的系数$\nu<0$使轮廓可以加快缩小；如果最初轮廓在目标对象的内部，加权区域的系数$\nu>0$使轮廓可以加快扩大。

通过解答PDE方程，得到最小化能量函数$\xi$最陡下降法的演化方程如下：
\begin{equation}\label{eq:3-14}
    \frac{\partial \phi}{\partial t}=\mu[\Delta\phi-div(\frac{\nabla \phi}{|\nabla \phi|})]+\mu\delta(\phi)div(g\frac{\nabla \phi}{|\nabla \phi|})+\nu g\delta(\phi)
\end{equation}
其中，$\Delta$为拉普拉斯算子（Laplacian Operator）。作为内能的$P(\phi)$，注意到其梯度流为：
\begin{equation}\label{eq:3-15}
    \Delta\phi-div(\frac{\nabla \phi}{|\nabla \phi|})=div[(1-\frac{1}{|\nabla \phi|})\nabla \phi]
\end{equation}
其中，$(1-1/|\nabla \phi|)$为扩散率因子。如果$|\nabla \phi|>1$，扩散率为正，这一项的效果通常是扩散，即使得$\phi$更均匀，因此可以减小$\phi$的梯度值$|\nabla \phi|$；相反，如果$|\nabla \phi|<1$，该项的效果是反向扩散，从而增大梯度。所以，由于作为惩罚项内能$P(\phi)$的存在，在$\phi$的演化过程中将自动保持一个近似于SDF的函数。因此，重新初始化过程可以完全消除。

与传统水平集方程相比，该模型方法的优势在于：1）无需重新初始化，加速曲线演化，缩短时间；2）水平集函数可以被初始化为相比PDE更容易计算的函数；3）此模型可以通过简单的有限差分实现，而不是传统水平集中的复杂的逆向方法（Upwind Scheme）计算。

\subsection{数值方法}
在实际的计算过程中，定义狄拉克（Dirac）函数$\delta$如下：
\begin{equation}\label{eq:3-16}
    \delta(x)=\left\{
        \begin{aligned}
            \frac{1}{2\xi}(1+cos\frac{\pi x}{\xi}),&\quad|x|\leq\xi;\\
            0,&\quad|x|>\xi.
        \end{aligned}
        \right.
\end{equation}

因为惩罚项所导致的扩散项（Diffusion Term）的存在，所以像在传统水平集方法中的逆向方法（Upwind Scheme）就没有必要了。取而代之的是，所有的空间偏微分${\partial \phi}/{\partial x}$和${\partial \phi}/{\partial y}$将近似为中心差分（Central Difference），时间偏微分${\partial \phi}/{\partial t}$近似为前向差分（Forward Difference）。则公式（\ref{eq:3-14}）通过近似差分的公式表达如下：
\begin{equation}\label{eq:3-17}
    \frac{\phi_{i,j}^{k+1} - \phi_{i,j}^{k}}{\tau}=L(\phi_{i,j}^k)
\end{equation}
其中，$L$是对公式（\ref{eq:3-14}）中等号右边部分的空间差分方法（Spatial Difference Scheme）的近似。上述公式的迭代方程为：
\begin{equation}\label{eq:3-18}
    \phi_{i,j}^{k+1}=\phi_{i,j}^{k}+\tau L(\phi_{i,j}^{k})
\end{equation}
\subsection{时间步长的选择}
在该方法的水平集实际计算中，时间步长$t$的选择比传统水平集方法的时间步长选择的自由性更大，例如从$0.1\text{到}100$在试验中都能达到不错的效果。

针对于迭代方程公式（\ref{eq:3-17}）而言，$t$的范围区间也是有限制的，$t\cdot\mu<1/4$才能保证水平集演化的稳定性。

选择较长的时间步长$t$可以加速演化，但是，如果时间步长$t$选择过大的时候，可能会导致在边界处出错（不精确）。所以，需要权衡两者：大的时间步长还是边界位置的精确性。一般而言，通常$t\leq10$适用于大多数图像。
\subsection{水平集函数的灵活初始化}
在传统水平集方法中，必须将水平集函数$\phi$初始化为一个符号距离函数$\phi_0$。如果初始化的水平集函数与SDF相差很大的时候，那么重新初始化过程也不可能将其重新初始化为一个SDF。在本文的方法中，不仅重新初始化过程可以完全消除，而且水平集函数$\phi$不再需要初始化为一个符号距离函数。

本章算法将提出以下函数作为初始化函数$\phi_0$。令$\Omega_0$为图像域$\Omega$的一个子集，$\partial \Omega$为在区域$\Omega_0$边界上的所有的点，即可以通过简单的形态学算子表示。这样，初始化函数$\phi_0$就可以定义为：
\begin{equation}\label{eq:3-19}
    \phi_0(x,y)=\left\{
        \begin{aligned}
            -d,&\quad (x,y) \in \Omega_0 - \partial \Omega_0;\\
            0,&\quad (x,y) \in \partial \Omega_0;\\
            d,&\quad \Omega - \Omega_0.
        \end{aligned}
    \right.
\end{equation}
其中，$d>0$为一个常数。且一般情况下，设定$d>2\xi$，$\xi$为正则化公式（\ref{eq:3-16}）的$\delta$函数中的宽度系数。

与由边界区域计算得到的SDF不同，上述的初始化水平集函数由在图像区域$\Omega$中的任意区域$\Omega_0$计算得到。这种基于区域的初始化水平集函数不仅计算效率高，而且在某些情况下还可以有灵活的应用。例如，如果感兴趣区域（ROI）可以以某种方式自动粗略地获得，如用阈值法，该方法就可以使用这个粗略获得区域$\Omega_0$构建初始化水平集函数$\Omega_0$。然后，初始水平集函数根据演化方程将稳定演化，其零水平集曲线正好是ROI的确切边界。

\section{实验结果}
本小节将对人工合成的图像和真实的超声图像进行分割实验，对比测地活动轮廓模型和改进后的活动轮廓之间的时间效率和分割结果的准确性。

首先，在合成图像的分割实验中，参数设定为：$\alpha=0.5, \beta = 0.2$，空间步长$h = 1$，时间步长$\Delta t = 0.1$。分割结果如图（\ref{fig:3-2}）所示。

\begin{figure}[htpb]

\centering
    \begin{minipage}[t]{0.5\textwidth}
    \label{fig:3-2(a)}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/GAC-hand-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    \label{fig:3-2(b)}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/GAC-hand-2600.jpg}\\(b)
    \end{minipage}
    \caption{(a)合成图像及初始化轮廓线，(b)合成图像分割结果}
\end{figure}
从上面人工合成图像（\ref{fig:3-2(a)}）和（\ref{fig:3-2(b)}）的分割结果可以看出：初始化轮廓线与图像真实边界相交的情况下，GAC模型也能得到良好的分割结果。GAC模型能得到这样效果的原因是，相对与几何活动轮廓模型增加的吸引因子项$\nabla g \cdot \nabla \phi$能够克服常量速度的作用，将目标边界内部的轮廓线吸引到真正的目标边界上来。GAC模型所具有的这样的优点，可以使得该模型应用于超声图像分割中有一定的优势，同时也存在一定的劣势，即在弱边界（或者噪声边界）上可能会使得该模型停止演化而达不到真正的边界。

同时，由于该模型在演化过程中需要重新初始化过程，因而该模型的在分割过程中所耗费的时间是比较大的。在接下来的真实超声图像分割实验中，将对比GAC模型和改进后的活动轮廓模型之间的分割效果以及时间效率。

在接下来的试验中，GAC模型的参数设置为和$\alpha=0.4$，空间步长$h = 1$，时间步长$\Delta t = 0.1$。并且在每$10$次迭代后进行重新初始化过程，以保证水平集函数符合SDF的特性。GAC模型改进后的模型的参数设置为$\lambda=6, \mu = 0.1, \nu = -3.0$，空间步长$h = 1$，时间步长$\tau = 2$。

\begin{figure}[htpb]
\centering
    \begin{minipage}[t]{0.5\textwidth}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/GAC-01-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/GAC-01-2200.jpg}\\(b)
    \end{minipage}
    \caption{(a)原始图像及初始化轮廓线，(b)GAC模型分割结果（迭代次数为2200，时间为68秒）}
\end{figure}

\begin{figure}[htpb]
\centering
    \begin{minipage}[t]{0.5\textwidth}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/LI-01-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/LI-01-900.jpg}\\(b)
    \end{minipage}
    \caption{(a)原始图像及初始化轮廓线，(b)改进模型分割结果（迭代次数为900，时间为12秒）}
\end{figure}

\begin{figure}[htpb]
\centering
    \begin{minipage}[t]{0.5\textwidth}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/GAC-08-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/GAC-08-1400.jpg}\\(b)
    \end{minipage}
    \caption{(a)原始图像及初始化轮廓线，(b)GAC模型分割结果（迭代次数为1400，时间为37秒）}

\end{figure}

\begin{figure}[htpb]
\centering
    \begin{minipage}[t]{0.5\textwidth}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/LI-08-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    \centering
    \includegraphics[scale=0.2]{figures/ch3/LI-08-400.jpg}\\(b)
    \end{minipage}
    \caption{(a)原始图像及初始化轮廓线，(b)改进模型分割结果（迭代次数为400，时间为5秒）}
\end{figure}

从图（\ref{fg:3-3}）和（\ref{fg:3-4}）可以看出，在相同的初始化轮廓线的条件下，改进后的模型比GAC模型不仅在时间上减小为GAC模型的六分之一，而且在分割效果上也比原始的模型好的多。图3-3和3-4为超声图像，图3-5和3-6是颈部血管的截面的超声图像。

从时间上考虑，GAC模型在图3-3和3-5中的迭代次数分别为2200次和1400次，因为传统的水平集方法中需要保证水平集函数的SDF特性，需要进行重新初始化过程。本实验中设定的重新初始化迭代次数为10，即每进行10次迭代后，将函数$\phi$进行重新初始化。然而，由于重新初始化过程是一个耗时的处理过程，所以在图像的分割中达到了68秒和37秒之多。而在改进后的模型中，达到图像中的分割结果所耗费时间约为GAC模型的六分之一，分别为12秒和5秒。

从分割结果的准确性考虑，如前面的分析，GAC模型在图像分割应用中有一定的优势，在超声图像的分割中也能达到一定的结果。然而由于超声图像本身的复杂性和分割的困难性，使得GAC在超声图像分割中并不能达到理想的效果。改进后的模型不仅保持了GAC模型中的吸引因子$\nabla g \cdot \nabla \phi$，还加入了面积项$A_g(\phi)$，将带分割区域的面积因素作为分割结果的一个判断因素。因为，重新初始化过程是对水平集函数的重新设置以满足SDF特性，会使得水平集函数偏离原来的真实水平集函数，所以会存在一定的误差。惩罚项的存在消去了重新初始化过程的必要性，不仅减少了程序所耗费的时间，而且也保证了水平集函数的准确性，所以能达到相对于GAC模型更准确的分割结果。

\section{本章小节}
本章主要讨论了基于边界信息的水平集活动轮廓模型，首先介绍由曲线演化和水平集方法相结合产生的第一个水平集活动轮廓模型：几何活动轮廓模型。由于几何活动轮廓模型只是曲线演化的一种应用（也称之为：轮廓演化模型），而在实际分割中并不广泛，其主要的意义就在于它将水平集方法应用于图像分割中，进而能自动地处理轮廓线的拓扑结构变化。

第二节通过分析Snake模型能量项推导出GAC模型，将能量函数与曲线演化相结合。本节解释了测地曲线与欧式曲线的联系和区别之处，还解释了GAC模型优于几何活动轮廓模型的原因，即边界吸引因子的存在。

针对包括GAC模型在内的传统水平集方法需要进行重新初始化过程的问题，第三节引入一种对GAC模型的改进模型，在模型中引入惩罚项和面积项（或加速项）。本节还讨论了该模型的数值实现以及时间步长$\tau$的选择范围等重要问题。

最后，通过实验验证了GAC模型的优点以及在超声图像分割中的应用，并比较了改进模型和GAC模型在时间和分割效果上的差别。本节通过对人工合成图像和真实超声图像的分割结果分析对比，改进后的模型在时间上缩短为GAC模型的六分之一，同时还提高了超声图像分割的准确性。

%--------------------------------------------第四章-------------------------
\chapter{基于区域信息活动轮廓模型的超声图像分割}
前一章所介绍的两种水平集活动轮廓模型有一个共同的特征，两种模型都利用待分割图像的图像边缘来分割图像，即两种模型都依赖图像的梯度$|\nabla u|$得到的边界函数$g$判断曲线演化的终止，都属于基于边界的图像分割模型\cite{chan2001active}。基于边界信息的水平集活动轮廓模型通常有以下缺点：

1.）对图像噪声比较敏感，通常在求图像边界之前需要对图像进行Gaussian滤波，然而当Gaussian滤波很效果很明显的时候，图像的边界信息也会变得模糊，而且基于边界信息的水平集活动轮廓模型仍会把强噪声误认为图像边缘，从而导致错误的分割；

2.）可能会发生边界泄漏现象，当基于边界信息的水平集活动轮廓模型的演化轮廓线与目标边界相交时候，通常会在常量速度的影响下离目标边界越来越远，导致最终的错误分割结果。虽然测地活动轮廓模型引入了吸引因子$\nabla g\cdot\nabla \phi$，在一定程度上避免了的边界泄漏现象的出现，但是如果目标边界上有相对比较弱的边缘时，测地活动轮廓模型也很难完全解决边界泄漏的问题；

3.）对初始轮廓线比较敏感，基于边界信息的水平集活动轮廓模型的初始化轮廓线需要被完全地设置在目标边界的内部或外部，否则初始化轮廓线与目标边界的相交部分会在常量速度的影响下远离真实的目标边界。

本章主要介绍另一类水平集活动轮廓模型，基于区域信息的水平集活动轮廓模型。该类模型利用图像的区域信息（主要是灰度信息，也可以是纹理信息或者统计信息等），可以克服基于边界信息的水平集活动轮廓模型所具有的缺点。Chan 和 Vese 首次提出了一种经典的基于区域的水平集活动轮廓模型——Chan-Vese 模型\cite{chan2001active}，该模型的停止条件是基于Mumford-Shah分割技术\cite{mumford1989optimal}的，因此该模型既可以检测有梯度的边界有可以检测无梯度的边界（Contours both with or without Gradient），比如光滑边界或者不连续边界。

本章结构安排如下：第一节介绍Chan-Vese模型以及该模型的数值解法，并将该方法应用与超声图像的分割中；第二节给出来Chan-Vese模型的一种改进模型，改进了Chan-Vese模型在不符合PC条件的图像分割效果；第三节提出了一种基于混合信息的活动轮廓模型，该模型综合利用边界信息和灰度信息对图像进行分割，并通过实验验证该方法的可行性和优点。

\section{C-V模型的原理与实现}
Chan-Vese模型由Mumford-Shah模型（M-S）演化而来，被认为是一种简化的M-S模型。Mumford-Shah模型\cite{mumford1989optimal}是一种变分模型，利用图像中目标边界曲线的特定规律实现图像分割。M-S模型假设图像由分片光滑区域构成，并用分片光滑函数表示灰度接近的同质区域，建立能量函数方程，并最小化该能量函数，将原始图像分割为（灰度）同质的连通区域。M-S模型的图像分割能量函数定义为：
\begin{equation}\label{eq:4-1}
    \begin{split}
        F^{MS}(u,C)& =\mu\cdot Length(C)\\
        & \quad +\lambda \int_\Omega{|u_0(x,y)-u(x,y)|} dxdy\\
        & \quad +\int_{\Omega\backslash C}{|\nabla u(x,y)|^2} dxdy
    \end{split}
\end{equation}
其中，$u_0$为给定图像，$\mu\text{和}\lambda$为正常数。通过最小化函数方程可以得到由曲线$C$包围的光滑区域$R_i$。一种M-S模型的简化形式是分片光滑的图像u，也就是说在每一个$\Omega\backslash C$区域内的连通区域$R_i$满足条件$u=\text{常数}c_i$。正如文献\cite{mumford1989optimal}指出在每个连通区域$R_i\text{上}c_i=average(u_0)$，这样的假设情况下，问题就简化为最小化分割问题（Minimal Partition Problem）。
\subsection{Chan-Vese模型}
Chan-Vese（C-V）模型\cite{chan2001active}是一种Mumford-Shah分割模型的简化形式。C-V模型假设图像$u$由两个近似分片光滑的均质区域构成，两个区域的灰度值分别为$u^i\text{和}u^o$，且待分割区域是灰度值为$u^i$的区域。根据此假设建立适当的能量项：
\begin{equation}
    \begin{split}
    F(C) &= F_1 (C)+F_2 (C)\\
    \quad &= \int_{inside C}{|u(x,y)-c_1 |^2} dxdy\\
    \quad &+ \int_{outside C}{|u(x,y)-c_2 |^2} dxdy
    \end{split}
\end{equation}
其中，$C$是任意可变曲线，常数$c_1\text{和}c_2$依赖于曲线$C$，分别为曲线内部和外部的灰度均值。不难发现，在这样的假设前提下，待分割目标边界$C_0$是最小化上述能量项的最终曲线，即：
\begin{equation}
    {\mathop{\mathrm{inf}}_{c_1, c_2,C}}\{F(c_1, c_2,C)\}\approx 0 \approx F_1(C+0)+F_2(C_0)
\end{equation}

考虑上述能量项在图（\ref{fig:4-1}）中各种情况下的能量，图（\ref{fig:4-1}）中A的情况下$F_1(C) > 0, F_2(C) \approx 0, F(C) > 0$，B情况下$F_1(C) \approx 0, F_2(C) > 0, F(C) > 0$，C情况下$F_1(C) > 0, F_2(C) > 0, F(C) > 0$，D情况下$F_1(C) \approx 0, F_2(C) \approx 0, F(C) = 0$。
\begin{figure}
    \centering
    \includegraphics[scale=0.5]{figures/ch4/CV-fitting.jpg}
    \caption{轮廓线与边界的各种位置情况}
    \label{fig:4-1}
\end{figure}



基于上述分析的能量项，再加上一些曲线的长度或者曲线包围区域的面积等规则化项，建立总的能量方程$F(c_1,c_2,C)$如下：
\begin{equation}
    \begin{split}
        F(c_1,c_2,C)& =\mu\cdot Length(C)+\nu\cdot Area(inside(C))\\
        & \quad +\lambda_1 \int_{inside C}{|u(x,y)-c_1|^2} dxdy\\
        & \quad +\lambda_2 \int_{inside C}{|u(x,y)-c_2|^2} dxdy
    \end{split}
\end{equation}
其中$\mu, \nu \geq 0; \lambda_1, \lambda_2 > 0$为常数参数。

最小化上述能量函数的过程可通过水平集方法实现。定义水平集函数$\phi$如下：
\begin{equation}
    \left\{
    \begin{aligned}
        C &=\partial \omega = \{(x,y)\in\Omega:\phi(x,y)=0\},\\
        inside(C) &= \omega = \{(x,y)\in\Omega:\phi(x,y)>0\},\\
        outside(C) &= \Omega \backslash \bar{\omega} = \{(x,y)\in\Omega:\phi(x,y)<0\}.
    \end{aligned}
    \right.
\end{equation}
通过引入如下Heaviside函数和Dirac函数将轮廓线$C$转化为用水平集函数$\phi$的方程表示为：
\begin{equation}
    \begin{split}
        F(c_1,c_2,\phi)
        & =\mu\int_\Omega{\delta(\phi(x,y))|\nabla \phi(x,y)|}dxdy\\
        & \quad +\nu\int_\Omega{H(\phi(x,y))}dxdy\\
        & \quad +\lambda_1\int_\Omega{(u(x,y)-c_1)^2H(\phi(x,y))}dxdy\\
        & \quad +\lambda_2\int_\Omega{(u(x,y)-c_2)^2(1-H(\phi(x,y)))}dxdy
    \end{split}
\end{equation}
最小化关于$c_1\text{和}c_2$的水平集能量函数$F(c_1,c_2,\phi)$且保持$\phi$不变，则$c_1\text{和}c_2$可表示为：
\begin{equation}
    \left\{
    \label{eq:4-7}
    \begin{aligned}
        c_1(\phi) &=\frac{\int_\Omega{u(x,y)H(\phi(x,y))}dxdy} {\int_\Omega{H(\phi(x,y))dxdy}}, \quad\text{如果}\int_\Omega{H(\phi(x,y))}>0;\\
        c_2(\phi) &=\frac{\int_\Omega{u(x,y)(1-H(\phi(x,y)))}dxdy} {\int_\Omega{(1-H(\phi(x,y)))dxdy}}, \quad\text{如果}\int_\Omega{(1-H(\phi(x,y)))>0}.
    \end{aligned}
    \right.
\end{equation}
相应的简化形式为：
\begin{equation}
    \left\{
    \label{eq:4-8}
    \begin{aligned}
        c_1(\phi) &= average(u) \text{如果}\{\phi \geq 0\};\\
        c_2(\phi) &= average(u) \text{如果}\{\phi < 0\}.
    \end{aligned}
    \right.
\end{equation}

用符号$C_0$表示该区域的边界。这样能得到，对于灰度值$u=u^i$的像素点就被$C_0$包括在内，灰度值$u=u^o$的就在$C_0$外部。

对于此变分活动轮廓模型的水平集方程，C-V模型将轮廓线$C$作为高维水平集的$\phi$的零水平集。对应于$\phi$的能量方程的演化方程为：
\begin{equation}\label{eq:4-9}
    \left\{
    \begin{split}
       \frac{\partial \phi}{\partial t}= \delta(\phi)[\mu div(\frac{\nabla \phi}{|\nabla \phi|}) - \nu - \lambda_1(u-c_1)^2\\
        \quad \quad+\lambda_2(u-c_2)^2] = 0,&\quad \text{在}(0, \infty)\times\Omega;\\
        \frac{\delta_\xi(\phi)}{|\nabla \phi|} \frac{\partial \phi}{\partial \vec n} = 0,&\quad \text{在边界}\partial\Omega\text{上};\\
        \phi(0,x,y)=\phi_0(x,y),&\quad \text{在}\Omega\text{内部}.
    \end{split}
    \right.
\end{equation}
其中，$\vec n$表示单位法线外向方向，$\partial \phi / \partial \vec n$表示$\phi$在边界处的法向微分。
\subsection{数值实现}
在实际的计算过程中Heaviside函数和Dirac函数选择有相对的自由性。如文献\cite{zhao1996variational}中引入的Heaviside函数，如下：
\begin{equation}
    H_{1,\xi}(z)=\left\{
        \begin{aligned}
            1,&\quad |z|>\xi;\\
            0,&\quad |z|<-\xi;\\
            \frac{1}{2}[1+\frac{z}{\xi}+\frac{1}{\pi}sin(\frac{\pi z}{\xi})],&\quad |z|\leq\xi.
        \end{aligned}
        \right.
\end{equation}

本文中使用文献\cite{chan2001active}中的正则化的H函数，如下：
\begin{equation}
    H_{2,\xi}=\frac{1}{2}(1+\frac{2}{\pi}arctan(\frac{z}{\xi}))
\end{equation}

\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.2]{figures/ch4/CV-H12.jpg}\\(a)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.2]{figures/ch4/CV-D12.jpg}\\(b)
%\label{5.4:side:b}
\end{minipage}
\caption{(a)不同的$H$函数$H_1, H_2$，(b)不同$H$函数下的$\delta$函数}
\end{figure}
两者的不同之处如图（\ref{fg:4-2}）所示。一个不同之处在于：$\delta_1$只在区间$-\xi, \xi$处不为$0$，而$\delta_2$在整个区间都不为$0$。因为能量为非凸包（允许有多个局部最小值，non-convex），所以，最终的结果可能依赖于初始化轮廓。使用$H_1$和$\delta_1$函数，此算法有时候计算一个局部的能量最小值，而使用$H_2$和$\delta_2$，此算法则更趋于计算一个全局的最小值。其中的一个原因是，在$H_1$条件下，关于$\phi$的欧拉-拉格朗日仅在$\phi=0$的零水平集的周围几个水平轮廓的局部范围内有效；而在$H_2$条件下，方程在每一个（级）水平集轮廓线上都有作用。如此，在实际应用中，可以得到一个与初始化轮廓线不相关的全局最小值；并且，也能够自动地探测到（图像）内部的轮廓。正如文献\cite{zhao1996variational}所述，为了能将演化扩展到$\phi$的所有的不同（级别的）水平集上，另一种可能的方法是将$\delta_0(\phi)$替换为$|\nabla\phi|$。本文依然使用$\delta_0(\phi)$，这样可以更接近于初始化最小问题（Initial Minimization Problem）。有关如何将演化扩展到$\phi$的所有的不同（级别的）水平集上，可以使用近似于$\delta_0$的$\delta_2$（$\delta_2$在任何地方都不为零）。

方程$\phi$离散化，将用到隐式有限差分方法（Finite Differences Implicit Scheme）。并使用以下标号：令$h$表示空间步长（Space Step），$\Delta t$表示时间步长（Time Step），对所有$1\leq -i, j\leq M, (x_i, y_j)=(ih,jh)$表示网格像素点。令$\phi_{i,j}^n=\phi(n\Delta t, x_i, y_j)$为$\phi(t,x,y)$的近似，$n\geq 0, \phi^0=\phi_0$，则有限差分为：
\begin{equation}
    \begin{aligned}
        \Delta_{\_}^x\phi_{i,j}& =\phi_{i,j}-\phi_{i-1,j},&
        \Delta_{+}^x\phi_{i,j}& =\phi_{i+1,j}-\phi_{i,j},\\
        \Delta_{\_}^y\phi_{i,j}& =\phi_{i,j}-\phi_{i,j-1},&
        \Delta_{+}^y\phi_{i,j}& =\phi_{i,j+1}-\phi_{i,j}.
    \end{aligned}
\end{equation}
根据上述定义，则算法过程为：对于已知的$\phi^n$，首先根据公式（\ref{eq:4-7}）分别计算$c_1$和$c_2$；然后，通过以下公式计算$\phi^{n+1}$：
\begin{equation}
    \begin{aligned}
        \frac{\phi_{i,j}^{n+1}-\phi_{i,j}^{n}}{\Delta t} &= \delta_h(\phi_{i,j}^{n}) \\
        \quad & [\frac{\mu}{h^2} \Delta_{\_}^x \cdot (\frac{\Delta_{+}\phi_{i,j}^n} {\sqrt{(\Delta_{+}^x\phi_{i,j}^n)^2/(h)^2 + (\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n})^2/(2h)^2}}) \\
        \quad &+\frac{\mu}{h^2} \Delta_{\_}^y \cdot (\frac{\Delta_{+}\phi_{i,j}^n} {\sqrt{(\Delta_{+}^y\phi_{i,j}^n)^2/(h)^2 + (\phi_{i,j+1}^{n}-\phi_{i,j-1}^{n})^2/(2h)^2}})\\
        \quad &-\nu-\lambda_1(u_{0,i,y}-c_1(\phi^n))^2 + \lambda_2(u_{0,i,y}-c_2(\phi^n))^2]
    \end{aligned}
\end{equation}
此线性化的方法是通过迭代方法实现的，详细部分可参考文献\cite{aubert1997variational}。

当使用到水平集和$\delta$函数时候，有一个标准的过程就是“重新初始化”，将$\phi$的零水平集曲线重新初始化为一个SDF函数，如文献\cite{zhao1996variational}所述，这样可以使得水平集函数不至于变得过于平滑（Flat），此过程也可以看做放缩或正则化过程。对于C-V模型方法而言，重新初始化过程不是必须的（Is Optional）；另一方面，此过程也不能过于强烈（Strong），因为，正如Fedkiw指出\cite{chan2001active}，这样可能会阻止内部轮廓增长。一般情况下，只对少数数值结果使用重新初始化过程，通过以下演化方程：
\begin{equation}
    \left \{
    \begin{aligned}
        \psi_\tau=sign(\phi(t))(1-|\nabla\psi|)\\
        \psi(0,\cdot)=\phi(t,\cdot)
    \end{aligned}
    \right.
\end{equation}
其中，$\phi(t)$为在时间$t$时候的$\phi$。下一个时刻的$\phi$是$\phi(t,\cdot)$，考虑各种情况，算法的主要流程可总结为如图（\ref{fig:4-3}）所示：
\begin{figure}[htpb]
    \thicklines
    \centering
    \includegraphics[scale=0.5]{figures/ch4/CV-pro.jpg}
    \caption{CV算法流程图}
    \label{fig:4-3}
\end{figure}

注意到，使用与时间相关的（Time-Dependent）关于$\phi$的PDE并不是至关重要的。从最小化问题得到的稳定问题（Stationary Problem）也可以使用相似的有限差分方法的数值解法解决。
\subsection{实验结果}
本节将通过三个实验测试C-V模型在超声图像和其他图像中的分割效果。首先，在符合C-V模型假设条件，即分片同质（PC）的条件下，且带分割目标为两个物体的图像中，设定参数为：$\mu=0.05*255^2, \lambda_1 = \lambda_2 = 1$，时间步长$\Delta t = 0.25$，分割结果如图（\ref{fig:4-6}）所示。
\begin{figure}[htpb]
\centering
    \begin{minipage}[t]{0.5\textwidth}
    %\label{fig:3-2(a)}
    \centering
    \includegraphics[scale=0.6]{figures/ch4/CV-09-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    %\label{fig:3-2(b)}
    \centering
    \includegraphics[scale=0.6]{figures/ch4/CV-09-240.jpg}\\(b)
    \end{minipage}
    \caption{(a)待分割图像及初始化轮廓线，(b)图像分割结果}
\end{figure}
从上面图像（\ref{fig:3-2(a)}）和（\ref{fig:3-2(b)}）的分割结果可以看出：CV模型的分割结果与初始化轮廓线的位置无关，只要待分割区域的灰度均一并且背景区域的灰度均一（PC条件），就能得到良好的分割结果，即能使得能量$F_1 (C)+F_2 (C)$达到最小的轮廓曲线。

将CV模型方法应用于超声图像的分割，且设定参数为：$\mu=0.0001*255^2, \lambda_1 = \lambda_2 = 1$，时间步长$\Delta t = 0.25$。分割结果如图（\ref{fig:4-6}）所示。

\begin{figure}[htpb]
\centering
    \begin{minipage}[t]{0.5\textwidth}
    %\label{fig:3-2(a)}
    \centering
    \includegraphics[scale=0.3]{figures/ch4/CV-us-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    %\label{fig:3-2(b)}
    \centering
    \includegraphics[scale=0.3]{figures/ch4/CV-us-400.jpg}\\(b)
    \end{minipage}
    \caption{(a)超声图像及初始化轮廓线，(b)超声图像分割结果}
\end{figure}

由于CV模型仅考虑全局的灰度分布，而且只有在满足模型的PC假设条件下才能达到较好的分割结果。由于超声图像的成像原理与其低分辨率、低对比度的特点，使得该方法在超声图像分割中的结果并不理想。如图（\ref{fig:4-5}）为前列腺超声图像，从图像（\ref{fig:4-5:b}）可看出，由于前列腺内部区域的不均一性，使得该方法将内部很多较亮的区域作为背景区域加以排除。

\begin{figure}[htpb]
\centering
    \begin{minipage}[t]{0.5\textwidth}
    %\label{fig:3-2(a)}
    \centering
    \includegraphics[scale=0.3]{figures/ch4/CV-dlg-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    %\label{fig:3-2(b)}
    \centering
    \includegraphics[scale=0.3]{figures/ch4/CV-dlg-400.jpg}\\(b)
    \end{minipage}
    \caption{(a)去斑点超声图像及初始化轮廓线，(b)去斑点超声图像分割结果}
\end{figure}
为了能近似满足CV模型的PC条件，首先将上述前列腺超声图像进行去斑点噪声处理，具体去斑点噪声可参考文献\cite{michailovich2006despeckling}。图（\ref{fig:4-6}）是对超声图像进行去斑点后的分割结果。从图中可以看出，在去斑点噪声后的图像中，相同的初始化轮廓线以及相同的迭代次数，在目标边界内部的亮点以及边界轮廓的毛刺都能被正确地划分，所达到的结果相对前者更加准确。

第三个实验是将C-V模型方法应用于数字减影血管造影（DSA）图像的分割中，如图fig:4-7所示：

\begin{figure}[htpb]
\centering
    \begin{minipage}[t]{0.5\textwidth}
    %\label{fig:3-2(a)}
    \centering
    \includegraphics[scale=0.6]{figures/ch4/CV-dsa1-000.jpg}\\(a)
    \end{minipage}%
    \begin{minipage}[t]{0.5\textwidth}
    %\label{fig:3-2(b)}
    \centering
    \includegraphics[scale=0.6]{figures/ch4/CV-dsa1-500.jpg}\\(b)
    \end{minipage}
    \caption{(a)血管图像及初始化轮廓线，(b)血管图像分割结果}
\end{figure}

图fig:4-7是数字减影血管造影图像的C-V模型分割结果。可以看出不仅图像中血管的边界模糊，而且图像中待分割血管的灰度也很不均一，在图像的左下角部分的血管的灰度与左上角的背景灰度比较接近。如果上述实验分析的原因相同，C-V模型很难能将全部的血管部分分割出来，而得不到正确的分割结果。

\section{C-V模型的改进模型RSF}
通过前一节的分析和实验说明，C-V模型适用范围较广并且也能得到良好的分割效果。同时C-V模型也存在着该模型本身的缺点，在其适用的图像中其演化速度慢；尤其是当待分割目标灰度不均一与背景灰度差较小的时候，该模型很难达到理想的分割效果。目前已有很多对C-V模型的改进算法，大致可以主要分为两类：快速算法\cite{shi2005fast,pan2006efficient,szpak2008further}和改进C-V模型在目标灰度不均一情况下的分割准确性\cite{sum2006novel,lin2008fast}。本节将给出一个基于局部区域灰度信息的活动轮廓模型，提高C-V模型在不符合PC条件下分割的准确性。

\subsection{RSF的实现}
为了能改进C-V模型加入局部灰度信息，需要引入一个包含以下特性的非负核函数$K$：

$1.) K(-u) = K(u);$

$2.) K(u)\geq K(v), \text{如果}|u|<|v|,\text{并且}lim_{|u|\rightarrow\infty}K(u)=0;$

$3.) \int k(x)dx = 1.$
其中，将第二个属性称为$K$核的“局部特性”。核函数$K$和它的局部特性将在此模型中起着重要的作用。

考虑一个给定的矢量（灰度或者彩色）图像，$d=1$为灰度图，$d=3$为彩色图像（RGB）。令$C$为一个给定图像$\Omega$中的闭合轮廓曲线，此轮廓线将$\Omega$分为两个区域$\Omega_1$和$\Omega_2$，分别表示轮廓线$C$的外部和内部区域。对一个给定的点$x\in\Omega$，定义“局部灰度拟合”项（Local Intensity Fitting Energy）：
\begin{equation}
    \xi_x^{Fit}(C, f_1(x), f_2(x)) = \sum\limits_{i=1}^2\lambda_i\int_{\Omega_i}{K(x-y)|I(y)-f_i(x)|^2}dy
\end{equation}
其中，$\lambda_i$为两个正数常量；$f(x)$分别为两个区域$\Omega_1$和$\Omega_2$的近似灰度均值；$I(y)$为以点$x$为中心的局部区域，其大小可以由核函数$K$来控制。因此，定义上述的LIF能量为点$x$在轮廓先$C$上的RSF能量。

核函数$K$的选择有较大的自由，只要它符合以上三个特性即可。因此，本文中设定核函数$K$为Gaussian核：
\begin{equation}
    K_\sigma(u)=\frac{1}{(2\pi)^{n/2}\sigma^n}e^{-|u|^2/2\sigma^2}
\end{equation}
其中$\sigma$为放缩因子。

接下来详细解释FIT（RSF）的意义：1. FIT是一个“加权均方误差”（Weighted Mean Square Error），是曲线$C$内部和外部灰度均值，并且在点$y$处的灰度值$I(y)$的权值为$K(x-y)$；2. 由于核函数局部特性的存在，（点$x$的）灰度值$I(y)$对FIT的贡献减小，并在远离中心点$x$处趋向于$0$。因此，能量FIT由在$x$邻域内的点$y$的灰度均值主导。特别地，当$y$远离$x$的时候，高斯核$K$迅速递减至$0$。粗略而言，当$|x-y|>3\sigma$，高斯核$K$有效性约为$0$。因此，仅在$y: |x-y|<3\sigma$的区域内的灰度才起作用。这样说来，这个拟合能量FIT是点$x$周围灰度的局部化。

在FIT中的适配能量是一个可放缩区域（Region-Scalable）：拟合值$f(x)$是以点$x$为中心的区域的灰度近似值，这个区域的大小可以由放缩因子$\sigma$控制。当放缩因子$\sigma$很小时，其所包含的区域就很小，相反，就会有更多的像素点被包含。注意到，如论文\cite{li2007implicit}所述，FIT充当的是局部拟合能量（Local Fitting Energy），其作用正好和C-V的（Piece-Constant，PC）模型中的全局拟合能量（Global Fitting Energy）相反。然而，称FIT为一个RSF更为合适，因为，FIT的灰度值不仅仅局限于一个小的局部区域。实际上，拟合能量FIT的灰度区域可以是任何大小的，可以是一个很小的邻域也可以是整个图像。区域放缩因子是唯一的也是新模型所需要的特征。给定一个中心点$x$，当轮廓线$C$恰为目标边界，并且拟合值$f(x)$最佳近似轮廓线$C$两边的灰度时，拟合能量FIT可以达到最小化。

为了得到整个目标边界，必须找到一个轮廓$C$，尽可能地减小在图像区域$\Omega$内对所有点$x$而言的拟合能量FIT。可以通过在图像区域，最大限度地降低所有中心点$x$的拟合能量FIT的积分达到，也就是最小化$\int{\xi_x^{Fit}(C,f_1 (x),f_2 (x))}dx$。此外，和在其他大部分的活动轮廓模型中一样，通过加入（或者，惩罚Penalize）轮廓线长度$C$，也有必要使得轮廓线$C$光滑。定义如下能量方程：
\begin{equation}
    \xi(C,f_1 (x),f_2 (x))=\int{\xi_x^{Fit}(C,f_1 (x),f_2 (x))}dx + \nu|C|
\end{equation}
这个能量函数根据轮廓线$C$定义，为了能自动处理拓扑变换，将采用水平集方法表示求解。

文献\cite{osher1988fronts}中，轮廓线$C$通过水平集函数（或者Lipschitz Function）$\phi$的零水平集表示。本文采用水平集函数$\phi$，令轮廓线$C$的外部和内部分别为正值和负值；$H$为传递函数（Heavisde Function），则能量函数$\int{\xi_x^{Fit}(C,f_1 (x),f_2 (x))}$可表示为：
\begin{equation}
\label{eq:4-18}
    \xi_x^{Fit}(C,f_1 (x),f_2 (x))=\sum\limits_{i=1}^2\lambda_i\int_{\Omega_i}{K(x-y)|I(y)-f_i(x)|^2}M_i(\phi(y)) dy
\end{equation}

其中，$M_1=H, M_2 = 1 - H$。因此，公式（\ref{eq:4-18}）可以重写为：
\begin{equation}
    \label{eq:4-19}
\begin{aligned}
    \xi(\phi, f_1, f_2) &= \sum\limits_{i=1}^2\lambda_i\int(\int{K(x-y)|I(y)-f_i(x)|^2}M_i(\phi(y)))dy dx \\
    \quad &+ \nu\int|\nabla H(\phi(x))|dx
\end{aligned}
\end{equation}
其中，最后一项$|\nabla H(\phi(x))|dx$计算$\phi$的零水平集的长度。

为了防止规则化（或，重新初始化）水平集函数$\phi$，本模型将引入一个水平集规则化项（Level Set Regularization Term）。如第三章模型一样，定义水平集规则化项为：
\begin{equation}\
    P(\phi)=\int {\frac{1}{2} {(|\nabla \phi|-1}^2)} dxdy
\end{equation}
如前所述，此项使得函数$\phi$在演化过程中更趋近于SDF。综合上述得到以下能量方程，最小化该能量的梯度流为新模型的水平集演化方程：
\begin{equation}
    F(\phi, f_1, f_2) = \xi_\xi(\phi, f_1, f_2) + \mu P(\phi)
\end{equation}

\subsection{能量最小化}
使用标准梯度下降（或最速下降）方法,最小化上述能量方程。对于一个固定的水平集函数$\phi$，最小化关于$f(x)$和$f(x)$的方$F(\phi,f_1,f_2 )$程$F(\phi,f_1,f_2 )$。通过计算方差（Variation），可以看出最小化能量方程$F(\phi,f_1,f_2 )$的$f1(x)$和$f2(x)$满足以下欧拉-拉格朗日（Euler-Lagrange）方程：
\begin{equation}
    \int K_\sigma(x-y)M_i^\epsilon(\phi(y))(I(y)-f_i(x))dy = 0, i=1,2
\end{equation}
得到：
\begin{equation}
    f_i(x)=\frac{K_\sigma(x)*[M_i^\epsilon(\phi(x))I(x)]}{K_\sigma(x)*M_i^\epsilon(\phi(x))}, i=1,2
\end{equation}
对于固定的$\phi$值，此$f_{1,2}(x)$能最小化能量方程$F(\phi,f_1,f_2)$。函数$f_{1,2}(x)$是在领域$x$内的加权灰度均值，其大小有放缩因子$\sigma$控制。注意到，公式中的分母始终为正数，因为$M_1=H>0$并且$M_2=1-H>0$。

保持$f_1$和$f_2$固定，最小化关于$\phi$的能量方程$F(\phi,f_1,f_2)$，通过使用标准梯度下降法得到的一下梯度流方程（Gradient Flow Equation）：
\begin{equation}
\label{eq:4-24}
\begin{aligned}
    \frac{\partial \phi}{\partial t} &= -\delta_\epsilon(\phi)(\lambda_1 e_1-\lambda_2 e_2) + \nu\delta_\epsilon(\phi)div(\frac{\nabla \phi}{|\nabla \phi|})\\
    \quad &- \mu(\nabla^2\phi-div(\frac{\nabla \phi}{|\nabla \phi|}))
\end{aligned}
\end{equation}
其中，$\delta_\epsilon$为平滑后的$\delta$，且$e_i$定义如下：
\begin{equation}
    e_i(x)=\int K_\sigma(y-x)|I(x)-f_i(y)|^2 dy, i=1,2
\end{equation}

公式（\ref{eq:4-22}）是新模型的水平集演化方程。第一项$-\delta_\epsilon(\phi)(\lambda_1 e_1-\lambda_2 e_2)$由数据拟合能量（Data Fitting Energy）推导出来，因此称为数据拟合项（Data Fitting Term）。这一项在新模型中的作用是：驱动活动轮廓模型趋向于目标边界。第二项$-\delta_\epsilon div(\nabla \phi / |\nabla \phi|$起到最短化轮廓长度（或者平滑零水平集曲线）的作用，它在维持轮廓规则化（Regularity）方面是必须的。所以，这一项被称为弧长项。第三项$-\mu(\nabla^2\phi-div(\nabla \phi / |\nabla \phi|)$称为水平集的规则化项，因为它在水平集演化过程中起到维持其规则（保持其为SDF）的作用。

\subsection{实验结果}
本节将通过实验说明RSF模型在不符合C-V模型的PC条件的图像分割中能达到比CV模型更好的结果。在本节实验中设定参数为：$\sigma = 3, \lambda_1 = \lambda_2 = 1$，时间步长$\Delta t = 0.1, \mu = 1, \nu = 0.001*255^2$。注意到，参数设定中放缩因子$\sigma$的取值较小，可以减小卷积运算。

\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.6]{figures/ch4/4-0.jpg}\\(a)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.5]{figures/ch4/4.jpg}\\(b)
%\label{5.4:side:b}
\end{minipage}
\caption{(a)原始图像及初始化轮廓线，(b)超声图像分割结果}
\end{figure}

\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.5]{figures/ch4/5-0.jpg}\\(a)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.5]{figures/ch4/5.jpg}\\(b)
%\label{5.4:side:b}
\end{minipage}
\caption{(a)原始图像及初始化轮廓线，(b)超声图像分割结果}
\end{figure}

从上面两幅血管图像的分割效果中可以看出，在CV模型不适合的待分割区域为非均匀区域的情况下，在目标血管图像的较暗的区域，CV模型都不能将其很好地分割出来，因为CV模型仅仅考虑全局的灰度信息而没有能将局部灰度信息作为分割的因素来考虑。RSF模型通过定义局部灰度拟合能量项，将局部的灰度变化作为分割算法的特征标准加入总的能量方程中，因此可以得到良好的分割结果。

\section{混合信息的活动轮廓模型}
GAC模型和其改进模型有效利用了图像的梯度信息，对图像边缘有良好的局部化分割效果。但是，由于此模型仅考虑了图像的局部梯度信息而没有考虑图像的全局信息，使得该模型对于边界模糊和边界无意义图像的分割效果不理想。

C-V模型不依赖于梯度而是利用区域信息对整个图像进行划分。该模型的能量函数包含图像区域信息，适用范围比基于梯度的要广，且对初始轮廓线位置不敏感。然而由于C-V模型的假设前提为该图像是两个均质区域，且只利用了全局灰度信息，显然在灰度层次比较丰富的图像中很难得到理想的分割效果。

本文提出综合上述两类模型的变分水平集活动轮廓模型，与单纯PDE驱动的水平集方法相比，变分水平集方法通常更方便；并且能更自然地将其他信息，比如区域信息和先验形状信息结合到能量方程中，并直接在水平集中公式化，能得到比较鲁棒的结果。

\subsection{惩罚项的引入}
在传统几何活动轮廓模型的水平集方法数值实现过程中，需要保证该演化的水平集函数$\phi$符合一个符号距离函数（SDF），即需要满足的一个特性是$|\nabla \phi|=1$，如3.3节分析，所以演化期间需要重新初始化，然而重新初始化是一个很耗时的计算过程，而且何时进行重新初始化也是一个很难的问题。新模型中引入惩罚项，该惩罚项使得水平集方程趋近一个SDF，从而完全消除耗时的重新初始化过程的必要性。引入3.3节中的惩罚项：
\begin{equation}
    P(\phi)=\int {\frac{1}{2} {(|\nabla \phi|-1}^2)} dxdy
\end{equation}
其中，$\phi$为水平集函数。该惩罚项作为衡量一个函数$\phi$趋近于SDF特征的标准，在演化过程中，规则化函数$\phi$使得其符合SDF的特征。引入惩罚项的原因在于新模型可以保证水平集函数符合SDF特性，并且能完全消除重新初始化的必要性。

\subsection{混合模型的建立}
为了能利用图像的区域信息，新模型引入C-V模型中的拟合能量项（Fitting Term）：
\begin{equation}
    \begin{split}
    C_r(C) &= F_1 (C)+F_2 (C)\\
    \quad &= \int_{inside C}{|u(x,y)-c_1 |^2} dxdy\\
    \quad &+ \int_{outside C}{|u(x,y)-c_2 |^2} dxdy
    \end{split}
\end{equation}
其中，$C$为任意变化的曲线，常量$c_1$和$c_2$根据$C$的变化而不同，分别为曲线$C$内部和外部的灰度均值。很明显，在这样的假设条件下，待分割目标的边界$C_0$为最小化能量项的最终曲线。

与C-V模型类似，新模型也是通过最小化拟合能量项达到分割的目的的。通过在新模型中加入一些规则化项，包括加权的曲线$C$的长度和曲线$C$内部的区域面积。本文提出能量函数方程$F(c_1, c_2, C)$如下：
\begin{equation}
    \begin{split}
        F(c_1,c_2,C)& =\lambda\cdot g\cdot Length(C)+\alpha \cdot g\cdot Area(inside(C))\\
        & \quad + \beta C_r(C)
    \end{split}
\end{equation}
其中，$g=1/(1+(\nabla G_\sigma *u)^2))$，是边界指示因子,参数$\lambda \geq 0$。公式中的第一项为加权曲线长度，其作用和GAC模型中的测地计算相似，将图像的边界$g$作为加权系数计算边界轮廓线的加权长度，显然只有当轮廓线位于目标边界时，此项才能到达最小值。第二项为加权面积项，有加速曲线演化的作用，当$\alpha$取正数时膨胀，相反就使得曲线逐渐收缩。

为了将上述关于曲线$C$的混合活动轮廓模型表示为水平集方程，通过未知变量$\phi$代替未知变量$C$。如此，能量方程$F(c_1, c_2, \phi)$可被写为：
\begin{equation}
    \begin{split}
        F(c_1,c_2,\phi)
        & =\lambda L_g(\phi) + \alpha C_g(\phi) + \beta C_r(\phi)\\
        & =\lambda\int_\Omega{\delta(\phi(x,y))|\nabla \phi(x,y)|}dxdy\\
        & \quad +\alpha\int_\Omega{H(\phi(x,y))}dxdy\\
        & \quad +\beta\int_\Omega{(u(x,y)-c_1)^2H(\phi(x,y))}dxdy\\
        & \quad +\beta\int_\Omega{(u(x,y)-c_2)^2(1-H(\phi(x,y)))}dxdy
    \end{split}
\end{equation}
同样的，新模型中的$c_1, c_2$的计算方法遵循C-V模型中的计算方法，具体公式为（\ref{eq:4-7}）或者（\ref{eq:4-8}）。

通过加入惩罚项，新模型的总的能量方程为：
\begin{equation}
    E(\phi)=\mu P(\phi) + \lambda L_g(\phi) + \alpha C_g(\phi) + \beta C_r(\phi)
\end{equation}
其中，$L_g$使得水平集函数的零水平集趋向于待分割目标的边界；$C_r$使得水平集函数的零水平集趋向于包含同质灰度区域；$C_g$加速曲线的演化过程。

通过解答PDE方程，得到最小化能量函数$E$最陡下降法的演化方程如下：
\begin{equation}
\begin{split}
    \frac{\partial \phi}{\partial t}&=\mu[\Delta\phi-div(\frac{\nabla \phi}{|\nabla \phi|})]\\
    \quad &+\lambda\delta(\phi)div(g\frac{\nabla \phi}{|\nabla \phi|})+\alpha g\delta(\phi)\\
    \quad &+ \beta \delta(\phi)((u(x,y)-c_1)^2-(u(x,y)-c_2)^2)
\end{split}
\end{equation}
由于该模型综合梯度和区域信息，曲线在模糊边界演化时使得轮廓线继续演化到合理的边界，因此能在弱边界和层次丰富的图像中得到良好的分割效果。

\subsection{实验结果与分析}
本节将通过实验说明基于混合信息的活动轮廓模型比单纯的基于边界信息或者基于区域信息的活动轮廓模型的分割效果有一定的优势。在本节实验中设定参数为：$\sigma = 3, \lambda_1 = \lambda_2 = 1$，时间步长$\Delta t = 0.1, \mu = 1, \nu = 0.001*255^2$。

本文通过一系列实验并对比各种模型对临床超声图像的分割效果。在以下各图中，(a)为初始化边界，(b)为基于边界的模型的分割效果，(c)为C-V模型的分割效果，(d)为本文提出的混合模型的分割效果。

\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-1a.jpg}\\(a)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-1b.jpg}\\(b)
%\label{5.4:side:b}
\end{minipage}
\caption{(a)原始图像及初始化轮廓线，(b)边界模型分割结果}
\end{figure}

\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.5]{figures/ch4/COM-1c.jpg}\\(c)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-1d.jpg}\\(d)
%\label{5.4:side:b}
\end{minipage}
\caption{(c)原始图像及初始化轮廓线，(d)混合模型分割结果}
\end{figure}


\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-2a.jpg}\\(a)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-2b.jpg}\\(b)
%\label{5.4:side:b}
\end{minipage}
\caption{(a)原始图像及初始化轮廓线，(b)边界模型分割结果}
\end{figure}

\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.4]{figures/ch4/COM-2c.jpg}\\(c)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-2d.jpg}\\(d)
%\label{5.4:side:b}
\end{minipage}
\caption{(c)原始图像及初始化轮廓线，(d)混合模型分割结果}
\end{figure}

图1和图2分别为对两组图像进行分割的结果。并且设置初始化轮廓在待分割区域的内部。从这两组图的对比中发现：(b)图的分割边界都相对比较光滑，边界的细节都未能很好地分割，(d)图比(b)图在边缘的细节方面有较大的提高；(d)图比(c)图在保证曲线的局部性方面也有较大的提高。


\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-3a.jpg}\\(a)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-3b.jpg}\\(b)
%\label{5.4:side:b}
\end{minipage}
\caption{(a)原始图像及初始化轮廓线，(b)边界模型分割结果}
\end{figure}

\begin{figure}[htpb]
\centering
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=0.5]{figures/ch4/COM-3c.jpg}\\(c)
%\label{5.4:side:a}
\end{minipage}%
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[scale=1]{figures/ch4/COM-3d.jpg}\\(d)
%\label{5.4:side:b}
\end{minipage}
\caption{(c)原始图像及初始化轮廓线，(d)混合模型分割结果}
\end{figure}

图1和图3分别为对相同图像的不同初始化轮廓线的分割结果。从这两组图的对比中发现：(d)图与(b)图相比的分割效果在不同的初始化轮廓线下能保持C-V模型所具有的特性，即对初始化边界不敏感，提高了分割效果的鲁棒性。


\section{本章小节}
本章主要讨论了基于区域信息的水平集活动轮廓模型，首先介绍了M-S分割模型理论，并由此引出首个基于区域的活动轮廓模型C-V模型，C-V模型是M-S模型的一种退化简单形式。从第一节的实验结果可以看出，该模型对初始化轮廓线不敏感，所以在实际应用比较广泛。然而，由于C-V模型有一个假设条件，即要求待分割图像符合PC条件，所以该模型在前景区域为非均一灰度的情况下很难达到理想的分割效果。

第二节针对C-V模型的PC条件假设前提的缺点，介绍了对该模型的一种改进模型RSF模型。RSF模型通过定义局部灰度拟合能量项，并且将局部的灰度变化作为分割算法的特征考虑，加入到总的能量方程中。实验表明，该模型在不符合PC假设条件的图像分割中能得到较之C-V模型更好的分割结果。

最后，本文提出了一种基于混合信息的活动轮廓模型图像分割方法。该方法在充分利用图像梯度信息的同时保留了区域信息来共同引导曲线的演化，因而相比单独依赖梯度或区域信息的方法在弱边界和层次丰富的超声图像能有较好的分割效果。实验表明该方法在超声图像分割上相对于基于边界的方法和基于区域C-V模型方法更为准确且具有鲁棒性。



%-----------------------------------------第五章-------------------------
\chapter{总结与展望}
医学影像在医学辅助诊断和检测发挥着重要的作用，借助于为医生提供人眼所不能直接观察到的疾病信息，能使得医生对疾病做出准确的判断。超声成像以其无创、无痛和无放射等特点越来越成为临床广泛应用的首选方式。对超声图像的精确分割或定位是HIFU治疗中的非常重要的环节，是后期重建、配准以及治疗能达到良好效果的前提。目前，国内外对超声图像分割已经有很多研究人员进行深入的研究，然而在基于活动轮廓模型的超声方面还未能系统全面。本文就这一主题进行了大量的研究工作，并将初步研究成果应用于实际的超声图像分割中，取得了一定的结果。

\section{本文主要研究工作}
讨论了基于边界信息的水平集活动轮廓模型，包括几何活动轮廓模型和测地活动轮廓模型等。首先简要介绍了几何模型的优点和缺点，并通过分析Snake模型能量项推导出GAC模型，将能量函数与曲线演化相结合。引入了测地曲线计算的概念，解释了测地曲线长度和欧式曲线长度的联系和区别，跟据两个模型的最终演化方程，解释了GAC模型优于几何活动轮廓模型的原因，即边界吸引因子的存在。

针对包括GAC模型在内的传统水平集方法需要进行重新初始化过程的问题，本文引入一种对GAC模型的改进模型，在模型中加入惩罚项和面积项（或加速项），从而完全消去了重新初始化的必要性，同时也加快了曲线的演化过程缩短了模型的演化时间。为了能保证结果的稳定性，讨论了该模型的数值实现以及时间步长的选择范围等重要问题。最后，通过实验验证了GAC模型的优点并将该方法应用于超声图像分割中，比较了改进模型和GAC模型在时间和分割效果上的差别。本文大量的实验验证了基于边界的水平集活动轮廓模型在超声分割中的可行性在超声分割中的可行性在超声分割中的可行性，同过经过对比分析实验结果，改进后的模型在时间上相对于GAC模型大大缩短，同时还提高了超声图像分割的准确性。

讨论了基于区域信息的水平集活动轮廓模型，首先简要介绍了M-S分割模型理论，并通过该模型解释了M-S模型的一种简化模型，首个基于区域的活动轮廓模型C-V模型。该模型对初始化轮廓线不敏感，实际应用比较广泛。然而，该模型在不符合C-V模型的PC假设条件的情况下，很难得到良好的分割结果。针对C-V模型的PC条件假设前提的缺点，介绍了对该模型的一种改进模型RSF模型。RSF模型通过定义局部灰度拟合能量项，并且将局部的灰度变化作为分割算法的特征考虑，加入到总的能量方程中，所以该模型在不符合PC假设条件的图像分割中能得到较之C-V模型更好的分割结果。

最后，本文提出了一种基于混合信息的活动轮廓模型图像分割方法。该方法在充分利用图像梯度信息的同时保留了区域信息来共同引导曲线的演化，因而相比单独依赖梯度或区域信息的方法在弱边界和层次丰富的超声图像能有较好的分割效果。实验表明该方法在超声图像分割上相对于基于边界的方法和基于区域C-V模型方法更为准确且具有鲁棒性。

\section{本文不足之处及研究方向}
本文所考虑的活动轮廓模型都是基于灰度信息和边界信息，但是，超声图像中与器官组织相关的纹理信息和形状信息，本文没有作为分割算法的特征考虑。

以后的研究工作可以从特定器官的纹理信息和活动轮廓模型相结合的方法来深入研究，例如在前列腺分割时，将纹理信息结合在活动轮廓模型。其他的改进模型可以结合形状模型或者灰度统计模型，进行深入的研究并应用于具体的超声图像分割中。

另外，本文中的初始化轮廓都是手工设定或者随机指定的，在初始化轮廓的设定中有很大的随机性。如果能结合分割算法的预处理，对分割算法最终分割精度的提高或许会有一定的帮助。

\end{Main} % 结束正文

%----------------------------------------------感 谢-------------------------
\begin{Thanks}
感谢……
\end{Thanks}

\bibliography{seuthesis}

%\begin{Appendix}
%  \chapter{第一个附录}
%  ……

%  \chapter{第二个附录}
%  ……
%\end{Appendix}

\newpage
\printindex % 索引

\begin{Resume}
作者简介
\end{Resume}

\backcover % 封底
\end{document}
